# Fit extreme value distribution to 
# sample maximum or minimum data
# Here we use the annual peak flow data
#Install and load the library fExtremes
library(fExtremes)
library(MASS)

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


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,1))

gevpdf=dgev(sort(peakflow), xi=zgev$par.ests[1], mu=zgev$par.ests[2], beta=zgev$par.ests[3])
gumbpdf=dgev(sort(peakflow), xi=0, mu=zgumb$par.ests[1], beta=zgumb$par.ests[2])
lognorpdf= dlnorm(sort(peakflow), meanlog=mean(log(peakflow)), sdlog=sd(log(peakflow)))

hist(peakflow, xlab="Peak flow", ylab="PDF", freq=FALSE,ylim=range(gevpdf,gumbpdf,lognorpdf))
points(peakflow,rep(0,length(peakflow))

lines(sort(peakflow), dgev(sort(peakflow), xi=zgev$par.ests[1], mu=zgev$par.ests[2], beta=zgev$par.ests[3]), col="red")

lines(sort(peakflow), dgev(sort(peakflow), xi=0, mu=zgumb$par.ests[1],
beta=zgumb$par.ests[2]), col="blue")

#Lognormal
lines(sort(peakflow), dlnorm(sort(peakflow), meanlog=mean(log(peakflow)), sdlog=sd(log(peakflow))), col="green")
box()
#Logpearson III = Gamma PDF on the Log of the data..

zgamma = fitdistr(log(peakflow), "gamma")
lp3pdf=dgamma(log(sort(peakflow)), shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])

hist(log(peakflow), freq=FALSE, xlab="Log (Peak Flow)", ylab="PDF",ylim=range(lp3pdf))
points(log(peakflow),rep(0,length(peakflow))

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 plot
sampquant = sort(peakflow)
N = length(peakflow)
sampprob = c(1:N)/(N+1)

modelquant = qgev(sampprob, xi=zgev$par.ests[1], mu=zgev$par.ests[2], beta=zgev$par.ests[3])
plot(sampquant, modelquant,xlab="SAmple Quantiles", ylab="Model Quantiles")
lines(sampquant,sampquant)