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) # Gumbel aa_gum <- (2*b1-b0)/log(2) cc_gum <- b0-0.5772*aa_gum vec_x <- cc_gum-aa_gum*log(-log(vec_p)) rr <- cor(vec_x,vec_q) # ss1 <- exp(-(vec_q-cc_gum)/aa_gum) 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_Gum.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","Gumbel") s1 <- sprintf(" n=%g\n",nn) s2 <- sprintf(" r=%.3f\n",rr) s3 <- sprintf(" a=%.3f\n",aa_gum) s4 <- sprintf(" c=%.3f\n",cc_gum) s5 <- 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) dev.off()