library(rjags)
library(MASS)


#### Fit GAmma PDF 
#### Estimate posterior distribution of Gamma PDF parameters
#### Plot the histogram of the posterior samples
#### Plot the PDFs from the posterior parameters 
#### overlay the PDF with the MLE of the original data

#### Fit the PDF for September flow at Lees Ferry 

test=matrix(scan("http://civil.colorado.edu/~balajir/r-session-files/Leesferry-mon-data.txt"),ncol=13,byrow=T)

flows=test[,2:13]/10^6     #get it into MAF

y = flows[,9]    ## Sep flow

### Fit GAMMA PDF - MLE

zz=fitdistr(y,"gamma")
shapefit = zz$estimate[1]
ratefit = zz$estimate[2]

### plot the histogram and overlay the PDF for visual

xeval=seq(max(0,min(y)-sd(y)),max(y)+sd(y),length=100)

xdensitygamma=dgamma(xeval,shape=zz$estimate[1],scale=1/zz$estimate[2])

zz=hist(y,freq=FALSE,plot=FALSE)
hist(y,xlab="Monthly rainfall", ylab="", freq=FALSE,
main="",ylim=range(c(zz$density, xdensitygamma)))
lines(xeval,xdensitygamma,col="blue")


############### Bayesian Parameter Estimation
## lambda = rate
### scale = 1/rate
## rr = shape


N = length(y)
model_string <- "
model{
  for(i in 1:N) {
    y[i] ~ dgamma(rr,lambda)
  }
  lambda ~ dgamma(0.01,0.01)  #prior
  rr ~ dgamma(0.01,0.01)
}
"
 model1.spec<-textConnection(model_string)

model <- jags.model(model1.spec, data = list(y = y,N=N), n.chains = 3, n.adapt= 10000)
update(model, 10000); # Burnin for 10000 samples
mcmc_samples <- coda.samples(model, variable.names=c("rr","lambda"), n.iter=20000)

postsamp = mcmc_samples[[1]]   #pull out the posterior simulations of lambda and rr

#### Posterior distribution of parameters
dev.new()
par(mfcol=c(2,1))
hist(postsamp[,1], main="Rate")   #lambda  = rate
abline(v=ratefit,col="red")

hist(postsamp[,2],main="shape")   #rr = shape	
abline(v=shapefit,col="red")

#Fit Gamma, plot PDFs from posterior parameters as grey scale
## overlay PDF from MLE

par(mfrow=c(1,1))

hist(y,freq=FALSE)

for(i in 1:500){
lines(xeval, dgamma(xeval, shape=postsamp[i,2], scale=1/postsamp[i,1]),col="grey")
}
lines(xeval,xdensitygamma,col="red",lwd=5)

points(y,rep(0,length(y)))
box()