# *****************************
# Histogram and GEV pdf
# *****************************
vec_x<- seq(1,400,1)

# Sapporo
d <- read.table('inp_RF_S.txt',skip=1) 
vec_qs <- sort(d[,1])
nn<-length(vec_qs)
vec_j <- c(1:nn)
vec_q<-vec_qs
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_ps <- 1/aa*(1-kk*(vec_x-cc)/aa)^(1/kk-1)*exp(-(1-kk*(vec_x-cc)/aa)^(1/kk))

# Maebashi
d <- read.table('inp_RF_M.txt',skip=1) 
vec_qm <- sort(d[,1])
nn<-length(vec_qm)
vec_j <- c(1:nn)
vec_q<-vec_qm
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_pm <- 1/aa*(1-kk*(vec_x-cc)/aa)^(1/kk-1)*exp(-(1-kk*(vec_x-cc)/aa)^(1/kk))

# Kyoto
d <- read.table('inp_RF_K.txt',skip=1) 
vec_qk <- sort(d[,1])
nn<-length(vec_qk)
vec_j <- c(1:nn)
vec_q<-vec_qk
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_pk <- 1/aa*(1-kk*(vec_x-cc)/aa)^(1/kk-1)*exp(-(1-kk*(vec_x-cc)/aa)^(1/kk))

# Plot
postscript(file="fig_R_GEV_HIST.eps", width=10,family="Helvetica")
par(mfrow=c(3,1)) 
hist(vec_qs,xlim=c(0,400),ylim=c(0,0.02),main="Sapporo",xlab="Observed data (mm/day)",probability=TRUE,breaks=seq(0,500,10),col="cyan")
lines(vec_x,vec_ps,type="l",lty=1,lwd=2,col="red")
hist(vec_qm,xlim=c(0,400),ylim=c(0,0.02),main="Maebashi",xlab="Observed data (mm/day)",probability=TRUE,breaks=seq(0,500,10),col="cyan")
lines(vec_x,vec_pm,type="l",lty=1,lwd=2,col="red")
hist(vec_qk,xlim=c(0,400),ylim=c(0,0.02),main="Kyoto",xlab="Observed data (mm/day)",probability=TRUE,breaks=seq(0,500,10),col="cyan")
lines(vec_x,vec_pk,type="l",lty=1,lwd=2,col="red")