flows=matrix(scan("http://civil.colorado.edu/~balajir/CVEN5454/R-sessions/sess2/Potomac-ann-max.txt"), ncol=2, byrow=T)
Years=flows[,1]
peakflow = flows[,2]
peakflow = peakflow / 10^3  #Thcfs

library(fExtremes)
#library(Rlab)

# mu1      # the threshhold
mu1 = 4
x = peakflow   #vector of data

  gpdmod = gpdFit(x,u=mu1)        #Fit GPD
  gpdpar = slot(gpdmod,"fit")    #get the parameters


  x1 = x[x >= mu1]    #get data exceeding the threshold


xeval = seq(mu1, (max(x1)+sd(x1)), length = 100)    #create a grid


gpdeval = dgpd(xeval,xi=gpdpar$par.ests[1],beta=gpdpar$par.ests[2],mu=mu1)




N = length(peakflow)

#ncls = round(log(N,2)+1)	#Stuge's formula for number of bins
#calculate monthly mean, variance and skews

N = length(x1)
#hplot(x1, xlab="Peak flow", ylab="PDF",nclass=ncls, ylim=range(gpdeval))

hist(x1,xlab="Peak flow", ylab="PDF", ylim=range(gpdeval),freq=FALSE)
points(x1,rep(0,N))
lines(xeval,gpdeval,col="red")


#### Q-Q plot

# get the sample CDF
weibp = rank(sort(x1))/(length(x1)+1)

## get the theoretical quantiles
## suppose we wish to check for GEV PDF
gpdquant = qgpd(weibp, xi=gpdpar$par.ests[1],beta=gpdpar$par.ests[2],mu=mu1)

plot(gpdquant, sort(x1), xlab="Theoretical Quantiles", ylab="Sample Quantiles")
lines(sort(x1), sort(x1))