d <- read.table('inp_RF_M.txt',skip=1) vec_q <- sort(d[,1]) nn <- length(vec_q) alp <- 0.4 #Cunnane formula vec_p <- ppoints(nn,alp) #plot(vec_q,vec_p) vec_j <- c(1:nn) b0 <- mean(vec_q) b1 <- sum((vec_j-1)*vec_q)/nn/(nn-1) b2 <- sum((vec_j-1)*(vec_j-2)*vec_q)/nn/(nn-1)/(nn-2) dd <- (2*b1-b0)/(3*b2-b0)-log(2)/log(3) kk <- 7.8590*dd+2.9554*dd*dd aa <- (2*b1-b0)/(1-2^(-kk))/gamma(1+kk)*kk cc <- b0-aa/kk*(1-gamma(1+kk)) vec_x <- cc+aa/kk*(1-(-log(vec_p))^kk) rr <- cor(vec_x,vec_q) # ss1 <- (1-kk*(vec_q-cc)/aa)^(1/kk) ss2 <- -log(vec_p) ww <- sum((ss1-ss2)^2) s99 <- -log(0.99) s01 <- -log(0.01) SLSC<- sqrt(ww/nn)/abs(s99-s01) # postscript(file="fig_R_HF_qq_GEV.eps", width=10,family="Helvetica") par(ps=18) par(pty="s") plot(vec_q,vec_x,log="xy",xlim=c(20,500),ylim=c(20,500),xlab="Observed Q",ylab="Estimated Q",pch=1,cex=1.2) abline(0,1) s0 <- sprintf(" %s\n","GEV") s1 <- sprintf(" n=%g\n",nn) s2 <- sprintf(" r=%.3f\n",rr) s3 <- sprintf(" k=%.3f\n",kk) s4 <- sprintf(" a=%.3f\n",aa) s5 <- sprintf(" c=%.3f\n",cc) s6 <- sprintf("SLSC=%.3f\n",SLSC) mtext(text=s0,side=3,line=-3.0,adj=0) mtext(text=s1,side=3,line=-4.5,adj=0) mtext(text=s2,side=3,line=-6.0,adj=0) mtext(text=s3,side=3,line=-7.5,adj=0) mtext(text=s4,side=3,line=-9.0,adj=0) mtext(text=s5,side=3,line=-10.5,adj=0) mtext(text=s6,side=3,line=-12.0,adj=0) dev.off()