#With the same data set as for the method of moments example, we fit a gamma distribution using the mle method.

illrain<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/ILLRAIN.txt",header=F,na.strings="*")

attach(illrain)

rain<-c(V1,V2,V3,V4,V5)
rain<-rain[!is.na(rain)]
rain

hist(rain,prob=T)

#The mle calculation is as follows:


227*log(mean(rain))

sum(log(rain))

(sum(log(rain))-227*log(mean(rain)))/227


f<-function(x) {log(x)-1.469825-digamma(x)}

uniroot(f,interval=c(0.01,3))

uniroot(f,interval=c(0.01,3))[1]

root<-as.numeric(uniroot(f,interval=c(0.01,3))[1])

root

alphahat<-root

lambdahat<-alphahat/mean(rain)


g<-function(x,alphahat,lambdahat){((lambdahat^alphahat)/gamma(alphahat))*x^(alphahat-1)*exp(-lambdahat*x)}


curve(g(x,.44079,1.964374),add=T)

#There is little practical difference between the method of moments and the method of maximum likelihood
#in the fit for this data set.

#Type or run

library(help=stats)

#to see a list of built in functions for statistical applications  which includes the function "uniroot".

#-----------------------------------------------------------------------------------------

#Bootstrapping to examine the variability of the mle estimate

resamp<-function(s){rgd<-rgamma(227,shape=.4408,rate=1.965)
                    
             f<-function(x) {227*log(x)-227*log(mean(rgd)) + sum(log(rgd)) -227*digamma(x)}
                    
                     alphahat<-uniroot(f,interval=c(0.01,3))[1]
                      
                     alphahat<-as.numeric(alphahat)
                     
                    lambdahat<-alphahat/mean(rgd)

                    alphahat

                    #lambdahat
 }

hist(sapply(1:1000,resamp))

alphadata<-sapply(1:1000,resamp)
hist(alphadata)
sqrt((1/1000)*sum((alphadata-mean(alphadata))^2))

#We change the above function slightly for the lambda computation


resamp<-function(s){rgd<-rgamma(227,shape=.4408,rate=1.965)
                    
             f<-function(x) {227*log(x)-227*log(mean(rgd)) + sum(log(rgd)) -227*digamma(x)}
                    
                     alphahat<-uniroot(f,interval=c(0.01,3))[1]
                      
                     alphahat<-as.numeric(alphahat)
                     
                    lambdahat<-alphahat/mean(rgd)

                    #alphahat

                    lambdahat
 }

lambdadata<-sapply(1:1000,resamp)
hist(lambdadata)
sqrt((1/1000)*sum((lambdadata-mean(lambdadata))^2))
----------------------------------------------------------------------------
#For comparison we now repeat some of the earlier example on the method of moments.
#Bootstrapping:(method of moments)

#Recall for the method of moments

varhat<-mean(rain^2)-mean(rain)^2

lambdahat<-mean(rain)/varhat

alphahat<-lambdahat* mean(rain)

alphahat

lambdahat

#The function we will use to do our bootstrapping is:

resamp<-function(x){rgd<-rgamma(227,shape=.3779,rate=1.6842)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    lambdastar
 }

lambdamoment<-hist(sapply(1:1000,resamp))



data<-sapply(1:1000,resamp)
slambdahat<-sqrt((1/1000)*sum((data-mean(data))^2))
slambdahat

#Now for the alphahat calculation: (changing the function only in the last line)

resamp<-function(x){rgd<-rgamma(227,shape=.3779,rate=1.6842)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    alphastar
 }

alphamoment<-hist(sapply(1:1000,resamp))

data<-sapply(1:1000,resamp)
salphahat<-sqrt((1/1000)*sum((data-mean(data))^2))
salphahat

-------------------------------------------------------------------
#Notice that the method of maximum likelihood estimates  have significantly less
# variance than those produced 
#by the method of moments for this data.

#This may be more clearly seen by comparing the histograms:

#First comparing the approximate distributions for the lambda estimates:

op<-par(no.readonly=TRUE)

par(mfcol=c(2,1),pin=c(1.8,.9))


hist(lambdadata,xlim=c(0,4))

plot(lambdamoment,xlim=c(0,4))
par(op)

#And now comparing the approximate distributions for the alpha estimates:


par(mfcol=c(2,1),pin=c(1.8,.9))


hist(alphadata,xlim=c(0,1))

plot(alphamoment,xlim=c(0,1))
par(op)

#Approximate confidence intervals(as in example E pg. 285):

#Using the estimated parameters from the maximum likelihood method, we derive
#an approximate  90% confidence interval for lambda_0.

o<-order(lambdadata)

lambdadata[o][50]

lambdadata[o][950]

lambdahat<- 1.965

deltalow<-lambdadata[o][50]- lambdahat

deltahigh<-lambdadata[o][950]-lambdahat

c(lambdahat-deltahigh,lambdahat-deltalow)

#And similarly for alpha_0:

o<-order(alphadata)

alphadata[o][50]

alphadata[o][950]

alphahat<- .4408

deltalow<-alphadata[o][50]- alphahat

deltahigh<-alphadata[o][950]-alphahat

c(alphahat-deltahigh,alphahat-deltalow)