############################################################
## Name: Joshua M. Tebbs
## Date: 22 Apr 2008
## Purpose: Fit gamma model to hurricane data using R
############################################################
 
# Enter data
y<-c(31,2.82,3.98,4.02,9.5,4.5,11.4,10.71,6.31,4.95,5.64,5.51,13.4,9.72,
6.47,10.16,4.21,11.6,4.75,6.85,6.25,3.42,11.8,0.8,3.69,3.1,22.22,7.43,5,
4.58,4.46,8,3.73,3.5,6.2,0.67)
 
## Second sample (uncentred) moment; needed for MOM
m2<-(1/36)*sum(y**2)
 
# MOM estimates (see Example 8.15 notes)
alpha.mom<-(mean(y))**2/(m2-(mean(y))**2)
beta.mom<-(m2-(mean(y))**2)/mean(y)
 
# Sufficient statistics
t1<-sum(log(y))
t2<-sum(y)
 
# Negative loglikelihood function (to be minimised)
# x1 = alpha
# x2 = beta
loglike<-function(x){
    x1<-x[1]
    x2<-x[2]
    36*log(gamma(x1))+36*x1*log(x2)-t1*(x1-1)+t2/x2
    }
 
# Use "optim" function to maximise the loglikelihood function
mle<-optim(par=c(alpha.mom,beta.mom),fn=loglike)
mle
 

# look at the qq-plot to assess the fit of the gamma model

plot(qgamma(ppoints(y),mle$par[1],1/mle$par[2]),sort(y),pch=16,

    xlab="gamma percentiles",ylab="observed values")