# Fit extreme value distribution to 
# sample maximum or minimum data
# Here we use the annual peak flow data

flows=matrix(scan("Potomac-ann-max.txt"), ncol=2, byrow=T)
Years=flows[,1]
peakflow = flows[,2]

flows=matrix(scan("http://civil.colorado.edu/~balajir/CVEN5454/R-sessions/sess2/LeesFerry-annmax-flow.txt"), ncol=4, byrow=T)
peakflow=flows[,4]
Years = 1884:2011

#Install and load the library fExtremes
library(fExtremes)

gevmodel = gevFit(peakflow)		#fit a Generalized Extreme Value Distribution
zgev=slot(gevmodel, "fit")

gumbelmodel = gumbelFit(peakflow)	#Fit a Gumbel distribution
zgumb = slot(gumbelmodel, "fit")


#### Plot them with histograms..
par(mfrow=c(2,2))
library(Rlab)
xeval = seq(max((min(peakflow)-sd(peakflow)),0), max(peakflow)+sd(peakflow),length=100)

gevdens = dgev(xeval, xi=zgev$par.ests[1], mu=zgev$par.ests[2], beta=zgev$par.ests[3])
gumbdens = dgev(xeval, xi=0, mu=zgumb$par.ests[1], beta=zgumb$par.ests[2])
lndens =  dlnorm(xeval, meanlog=mean(log(peakflow)), sdlog=sd(log(peakflow)))
#hist(peakflow, xlab="Peak flow", ylab="PDF", probability=T)

ncls = log(length(peakflow),2)+1
hplot(peakflow, nclass=ncls,xlab="Peak flow", ylab="PDF", ylim=range(gevdens, gumbdens, lndens),xlim=range(xeval))
lines(xeval, gevdens, col="red")
lines(xeval,gumbdens, col="blue")
lines(xeval,lndens, col="green")

box()

#Logpearson III = Gamma PDF on the Log of the data..
library(MASS)
zgamma = fitdistr(log(peakflow), "gamma")

hplot(log(peakflow), nclass=ncls, xlab="Log (Peak Flow)", ylab="PDF")

lines(log(sort(peakflow)), dgamma(log(sort(peakflow)), shape=zgamma$estimate[1], scale=1/zgamma$estimate[2]),col="purple")

box()
################

#Find the Quantiles 50-year, 100-year and 500-year return period.
# for GEV
gevquant = qgev(c((1-1/50), (1-1/100), (1-1/500)), xi=zgev$par.ests[1], mu=zgev$par.ests[2], beta=zgev$par.ests[3])

### You can similarly compute the return periods for other distributions..

### Q-Q plots..
par(mfrow=c(2,2))

#modquant = qgev(ppoints(length(peakflow),a=0),xi=zgev$par.ests[1],mu=zgev$par.ests[2], beta=zgev$par.ests[3])
modquant = qgev(ppoints(length(peakflow)),xi=zgev$par.ests[1],mu=zgev$par.ests[2], beta=zgev$par.ests[3])
plot(modquant, sort(peakflow),xlab="Theoretical Quantile", ylab="Sample Quantile",ylim=range(modquant,peakflow), xlim=range(modquant, peakflow))
lines(sort(peakflow),sort(peakflow),col="red")
title(main="Q-Q plot from GEV")

## Gumbel
modquant = qgev(ppoints(length(peakflow),a=0),xi=0,mu=zgumb$par.ests[1], beta=zgumb$par.ests[2])
plot(modquant, sort(peakflow),xlab="Theoretical Quantile", ylab="Sample Quantile",ylim=range(modquant,peakflow), xlim=range(modquant, peakflow))
lines(sort(peakflow),sort(peakflow),col="red")
title(main="Q-Q plot from Gumbel")

# Lognormal
modquant = qlnorm(ppoints(length(peakflow),a=0),meanlog=mean(log(peakflow)), sdlog=sd(log(peakflow)))
plot(modquant, sort(peakflow),xlab="Theoretical Quantile", ylab="Sample Quantile",ylim=range(modquant,peakflow), xlim=range(modquant, peakflow))
lines(sort(peakflow),sort(peakflow),col="red")
title(main="Q-Q plot from LogNormal")


# LogPearson type 3
modquant = qgamma(ppoints(length(peakflow),a=0),shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])
plot(modquant, log(sort(peakflow)),xlab="Theoretical Quantile", ylab="Sample Quantile",ylim=range(modquant,log(peakflow)), xlim=range(modquant, log(peakflow)))
lines(log(sort(peakflow)),log(sort(peakflow)),col="red")
title(main="Q-Q plot from LogPearson III")