/******Problem 46 ******/ A- Make a histogram of the 210 values of ti. Does it appear that a gamma distribution would be a plausible model to ???t?

Data has been import by import option.

#Reading the dataset
whale=read.table("C:/Users/chink/Google Drive/Computational Statistics/4th Quater/Mathematical Statistics I/Homework/Homework_03/whale.txt")
#Creating Breaks 
whale=whale$V1
br = seq(0,5,0.25)
hist(whale,breaks=br,main='whale data',probability = T,plot=T,col='green')

Yes, Histogram is look like gammma distribution.

out.mme = matrix(0,2,5)

out.mle = matrix(0,2,7)
  1. Fit the parameters of the gamma distribution by the method of moments.

More details about Gamma distribution http://thecoatlessprofessor.com/statistics/proofs-of-gamma-distribution/

#Finding parameters from gamma distribution
n=length(whale)         # sample size     
X.bar=mean(whale)        #Sample Mean
sigmasq.hat = ((n-1)/n)*var(whale)  # (1/n)sum(x_i - x.bar)^2
lambda.hat.mme = X.bar/sigmasq.hat
lambda.hat.mme
## [1] 1.329121
alpha.hat.mme = X.bar^2/sigmasq.hat
alpha.hat.mme
## [1] 0.8002543
#Keeping MME in matrix
out.mme[1,1] = alpha.hat.mme
out.mme[2,1] = lambda.hat.mme
  1. Fit the parameters of the gamma distribution by maximum likelihood. How do these values compare to those found before.
# solve the nonlinear equation
alpha.d = function(a){
    n*log(a) - n*(log(mean(whale))) + sum(log(whale)) - n*digamma(a)
}
a.min = 0.1
a.max = 5
alpha.hat.mle = uniroot(alpha.d, c(a.min,a.max))

alpha.hat.mle = alpha.hat.mle$root
alpha.hat.mle
## [1] 1.611355
lambda.hat.mle = alpha.hat.mle / mean(whale)
lambda.hat.mle
## [1] 2.676257
out.mle[1,1] = alpha.hat.mle
out.mle[2,1] = lambda.hat.mle
   Lambda.hat     Alpha.hat

MME 1.329121 0.8002543 MLE 2.676257 1.611355

  1. Plot the two gamma densities on top of the histogram. Do the ???ts look reasonable?
# plot the fitted density using mme
par(new=T)
## Warning in par(new = T): calling par(new=TRUE) with no plot
x = seq(0.01,5,0.01)
y = dgamma(x,alpha.hat.mme,lambda.hat.mme)
y.max = 1.8 
plot(x,y,type='l',ylim=c(0,y.max),col="green")

# plot the fitted density using mle
par(new=T)
## Warning in par(new = T): calling par(new=TRUE) with no plot
x = seq(0.01,5,0.01)
y = dgamma(x,alpha.hat.mle,lambda.hat.mle)
y.max = 1.8
plot(x,y,type='l',lty=4,ylim=c(0,y.max),col="green")

In the mme plot at x=0, y value is approching towards infinity while in the mle plot, at x=0, y value is approcting towards the0.

  1. Estimate the sampling distributions and the standard errors of the parameters ???t by the method of moments by using the bootstrap.
B = 1000

y = matrix(0,B,n)

lambda.star.mme = numeric(B)
alpha.star.mme = numeric(B)

for(i in 1:B){
    y[i,] = rgamma(n, alpha.hat.mme,lambda.hat.mme)
    x.bar = mean(y[i,])
    sigmasq.hat = ((n-1)/n)*var(y[i,]) 
    lambda.star.mme[i] = x.bar/sigmasq.hat
    alpha.star.mme[i] = x.bar^2/sigmasq.hat
}

br = seq(0,5,0.10)  # create breaks for the histogram
hist(alpha.star.mme,breaks=br,probability = T,plot=T,col="green") 

alpha.star.mean = mean(alpha.star.mme)
alpha.star.mean
## [1] 0.8249198
alpha.star.sd = sqrt(var(alpha.star.mme))
alpha.star.sd
## [1] 0.1190054
out.mme[1,2] = alpha.star.mean
out.mme[1,3] = alpha.star.sd

alpha.LCL = quantile(alpha.star.mme, 0.025)
alpha.LCL
##      2.5% 
## 0.6109046
alpha.UCL = quantile(alpha.star.mme, 0.975)
alpha.UCL
##    97.5% 
## 1.075058
out.mme[1,4] = alpha.LCL
out.mme[1,5] = alpha.UCL
hist(lambda.star.mme,breaks=br,probability = T,plot=T,col="green")

lambda.star.mean = mean(lambda.star.mme)
lambda.star.mean
## [1] 1.379045
lambda.star.sd = sqrt(var(lambda.star.mme))
lambda.star.sd
## [1] 0.2191114
out.mme[2,2] = lambda.star.mean
out.mme[2,3] = lambda.star.sd
lambda.LCL = quantile(lambda.star.mme, 0.025)
lambda.LCL
##      2.5% 
## 0.9584404
lambda.UCL = quantile(lambda.star.mme, 0.975)
lambda.UCL
##    97.5% 
## 1.835086
out.mme[2,4] = lambda.LCL
out.mme[2,5] = lambda.UCL
  1. Estimate the sampling distributions and the standard errors of the parameters ???t by maximum likelihood by using the bootstrap. How do they compare to the results found previously?
B = 1000

y = matrix(0,B,n)

lambda.star.mle = numeric(B)
alpha.star.mle = numeric(B)
a.min = 0.5
a.max = 5

for(i in 1:B){
    y[i,] = rgamma(n, alpha.hat.mle,lambda.hat.mle)
    whale = y[i,]
    alpha.star.mle[i] = uniroot(alpha.d, c(a.min,a.max))$root   # alpha
    lambda.star.mle[i] = alpha.star.mle[i] / mean(whale)        # lambda
}

br = seq(0,5,0.05)  # create breaks for the histogram
hist(alpha.star.mle,breaks=br,probability = T,plot=T,col="green") 

alpha.star.mean = mean(alpha.star.mle)
alpha.star.mean
## [1] 1.642144
alpha.star.sd = sqrt(var(alpha.star.mle))
alpha.star.sd
## [1] 0.15207
out.mle[1,2] = alpha.star.mean
out.mle[1,3] = alpha.star.sd
alpha.LCL = quantile(alpha.star.mle, 0.025)
alpha.LCL
##     2.5% 
## 1.374308
alpha.UCL = quantile(alpha.star.mle, 0.975)
alpha.UCL
##    97.5% 
## 1.972585
out.mle[1,4] = alpha.LCL
out.mle[1,5] = alpha.UCL
hist(lambda.star.mle,breaks=br,probability = T,plot=T,col="green")

lambda.star.mean = mean(lambda.star.mle)
lambda.star.mean
## [1] 2.730994
lambda.star.sd = sqrt(var(lambda.star.mle))
lambda.star.sd
## [1] 0.2925581
out.mle[2,2] = lambda.star.mean
out.mle[2,3] = lambda.star.sd
lambda.LCL = quantile(lambda.star.mle, 0.025)
lambda.LCL
##     2.5% 
## 2.202238
lambda.UCL = quantile(lambda.star.mle, 0.975)
lambda.UCL
##    97.5% 
## 3.337758
out.mle[2,4] = lambda.LCL
out.mle[2,5] = lambda.UCL
  1. Find approximate con???dence intervals for the parameters estimated by maximum likelihood.
# alpha 

delta.low = quantile(alpha.star.mle,0.025) - alpha.star.mean 
delta.up = quantile(alpha.star.mle,0.975) - alpha.star.mean

mle.LCL = alpha.hat.mle - delta.up 
mle.LCL
##    97.5% 
## 1.280915
mle.UCL = alpha.hat.mle - delta.low
mle.UCL
##     2.5% 
## 1.879191
out.mle[1,6] = mle.LCL
out.mle[1,7] = mle.UCL
# lambda 

delta.low = quantile(lambda.star.mle,0.025) - lambda.star.mean 
delta.up = quantile(lambda.star.mle,0.975) - lambda.star.mean

mle.LCL = lambda.hat.mle - delta.up 
mle.LCL
##    97.5% 
## 2.069494
mle.UCL = lambda.hat.mle - delta.low
mle.UCL
##     2.5% 
## 3.205013
out.mle[2,6] = mle.LCL
out.mle[2,7] = mle.UCL

Confidence intervals for the alpha.hat is ( 1.271732 , 1.880815 ). confidence interval for the lambda.hat is ( 2.005162 , 3.193238 )

dimnames(out.mme) = list(c("alpha","lambda"),c("MME", "Bootstrap.mean", "Bootstrap.sd", "quant.LCL", "quant.UCL"))
dimnames(out.mle) = list(c("alpha","lambda"),c("MLE", "Bootstrap.mean", "Bootstrap.sd", "quant.LCL", "quant.UCL", "approx.LCL", "approx.UCL"))

out.mme
##              MME Bootstrap.mean Bootstrap.sd quant.LCL quant.UCL
## alpha  0.8002543      0.8249198    0.1190054 0.6109046  1.075058
## lambda 1.3291210      1.3790454    0.2191114 0.9584404  1.835086
out.mle
##             MLE Bootstrap.mean Bootstrap.sd quant.LCL quant.UCL approx.LCL
## alpha  1.611355       1.642144    0.1520700  1.374308  1.972585   1.280915
## lambda 2.676257       2.730994    0.2925581  2.202238  3.337758   2.069494
##        approx.UCL
## alpha    1.879191
## lambda   3.205013