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)

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)

vec_x <- seq(vec_q[1],vec_q[nn],1)

# Gumbel
aa_gum <- (2*b1-b0)/log(2)
cc_gum <- b0-0.5772*aa_gum
vec_y_gum <- -log(-log(exp(-exp(-(vec_x-cc_gum)/aa_gum))))

# GEV
dd <- (2*b1-b0)/(3*b2-b0)-log(2)/log(3)
kk_gev <- 7.8590*dd+2.9554*dd*dd
aa_gev <- (2*b1-b0)/(1-2^(-kk_gev))/gamma(1+kk_gev)*kk_gev
cc_gev <- b0-aa_gev/kk_gev*(1-gamma(1+kk_gev))
vec_y_gev <- -log(-log(exp(-(1-kk_gev*(vec_x-cc_gev)/aa_gev)^(1/kk_gev))))

postscript(file="fig_R_GUM_paper.eps", width=10,family="Helvetica")
par(ps=16)
par(pty="s")
plot(vec_q,-log(-log(vec_p)),xaxt="n",yaxt="n",xlab="",ylab="",
xlim=c(0,500),ylim=c(-2,7),
main="Gumbel probability paper",pch=1,cex=1.2)
title(xlab="Observed Q", mgp=c(3,1,0))
title(ylab="Non-exceedance probability", mgp=c(4,1,0))

prob<-c(0.001,0.01,0.1,0.5,0.9,0.99,0.999)
yax<- -log(-log(prob))
xax<-seq(0,500,100)
lines(vec_x,vec_y_gum,type="l",lty=1,lwd=1,col="red")
lines(vec_x,vec_y_gev,type="l",lty=1,lwd=1,col="blue")
abline(h=yax,v=xax,col="gray")
axis(2,at=yax,labels=prob,cex.axis=1.0,las=2)
axis(1,at=xax,labels=xax,cex.axis=1.0)
par(bg="white")
legend("topleft",legend=c("Cunnane","Gumbel","GEV"),
col=c("black","red","blue"),lty=c(0,1,1),pch=c(1,-1,-1),
cex=0.9,pt.cex=1.2,y.intersp=1.5)
dev.off()