arrivals<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/gamma-arrivals.txt",header=F)

names(arrivals)

attach(arrivals)

hist(V1,prob=T)

A gamma distribution is a plausible model since the histogram 
is skewed to the left.

a<-log(mean(V1))-mean(log(V1))

f<-function(x) {log(x)-a-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

alphahatmle<-root

lambdahatmle<-alphahatmle/mean(V1)


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


curve(g(x,1.026,.01284),add=T) #m.l.e

curve(g(x,1.0124,.01266),add=T) #moments


Bootstrapping for standard error:

length(V1)

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

                    #alphahat

                    lambdahat
 }


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

#.02024 in one determination.

#We change the above function slightly for the lambda computation.
#Move the comment mark to the appropriate parameter before recompiling it.



lambdadata<-sapply(1:1000,resamp)
hist(lambdadata)
sqrt((1/1000)*sum((lambdadata-mean(lambdadata))^2))
#.000328  in one determination.



Method of Moments:

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

lambdahat<-mean(V1)/varhat

alphahat<-lambdahat* mean(V1)

alphahat

lambdahat

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

resamp<-function(x){rgd<-rgamma(length(V1),shape=alphahat,rate=lambdahat)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    lambdastar
 }




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

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

resamp<-function(x){rgd<-rgamma(length(V1),shape=alphahat,rate=lambdahat)
                    mean(rgd^2)-mean(rgd)^2->var
                    lambdastar<-mean(rgd)/var
                    alphastar<-lambdastar*mean(rgd)
                    alphastar
 }



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

Bootstrapping for confidence intervals:

m.l.e:

o<-order(lambdadata)

lambdadata[o][50]

lambdadata[o][950]


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

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

c(lambdahatmle-deltahigh,lambdahatmle-deltalow)

#And similarly for alpha_0:

o<-order(alphadata)

alphadata[o][50]

alphadata[o][950]


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

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

c(alphahatmle-deltahigh,alphahatmle-deltalow)

method of moments:

o<-order(lambdamomentdata)

lambdamomentdata[o][50]

lambdamomentdata[o][950]


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

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

c(lambdahat-deltahigh,lambdahat-deltalow)

#[1] 0.01193909 0.01343938

#And similarly for alpha_0:

o<-order(alphadata)

alphadata[o][50]

alphadata[o][950]


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

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

c(alphahat-deltahigh,alphahat-deltalow)

#[1] 0.9771597 1.0433721