################################################### gevpar=function(x){ n=length(x) sx=sort(x) l1=0; l2=0; l3=0; for (i in 1:n) { l1=l1+sx[i] l2=l2+sx[i]*(2*i-n-1) l3=l3+sx[i]*(6*i*(i-n-1)+n*n+3*n+2) } l1=l1/n; l2=l2/n/(n-1); l3=l3/n/(n-1)/(n-2) t3=l3/l2 zz=2/(3+t3)-log(2)/log(3) kk=7.859*zz+2.9554*(zz^2) aa=l2*kk/((1-2^(-kk))*gamma(kk+1)) ee=l1+aa*(gamma(kk+1)-1)/kk if ((kk<0)&(min(x)<=ee+aa/kk)) kk=aa/(min(x)-ee) if ((kk>0)&(max(x)>=ee+aa/kk)) kk=aa/(max(x)-ee) par=0 par[1]=ee par[2]=aa par[3]=kk gevpar=par } ########################################################## data=read.table("tmax_LOCID000035.txt") data=as.matrix(data) X=data[,2] nt=length(X) rp=30 prob.rp=1-1/rp ee=gevpar(X)[1]; aa=gevpar(X)[2]; kk=gevpar(X)[3] sX=sort(X) f.emp=1:nt/(nt+1) f.gev=exp( -(1-(sX-ee)*kk/aa)^(1/kk) ) rv=ee+aa*(1- (-log(1-1/rp))^kk)/kk plot(f.gev,sX,type="l",main="Paris, Tmax",ylim=c(29,41)) points(f.emp,sX) lines(c(prob.rp,prob.rp),c(min(X)-1,rv),col="red") lines(c(-0.1,prob.rp),c(rv,rv),col="red") axis(side=3,at=c(0.5,prob.rp),c(2,30)) data00=read.table("tmax_00_Paris.txt") data00=as.matrix(data00) X00=data00[,2] nt=length(X00) rp=30 prob.rp=1-1/rp ee00=gevpar(X00)[1]; aa00=gevpar(X00)[2]; kk00=gevpar(X00)[3] sX00=sort(X00) f.emp=1:nt/(nt+1) f.gev00=exp( -(1-(sX00-ee00)*kk00/aa00)^(1/kk00) ) rv00=ee00+aa00*(1- (-log(1-1/rp))^kk00)/kk00 data03=read.table("tmax_03_Paris.txt") data03=as.matrix(data03) X03=data03[,2] nt=length(X03) rp=30 prob.rp=1-1/rp ee03=gevpar(X03)[1]; aa03=gevpar(X03)[2]; kk03=gevpar(X03)[3] sX03=sort(X03) f.emp=1:nt/(nt+1) f.gev03=exp( -(1-(sX03-ee03)*kk03/aa03)^(1/kk03) ) prob03=exp( -(1-(rv00-ee03)*kk03/aa03)^(1/kk03) ) T03=1/(1-prob03) plot(f.gev00,sX00,type="l",main="Paris, Tmax",ylim=c(36,44)) lines(f.gev03,sX03,type="l",lty="dashed")