STA 506 Assignment 5: Maximum Likelihood Estimation (MLE)

1 Problem 1

For this problem, we are asked to derive the general log-likeihood function of the gamma distribution to estimate its parameters, \(\alpha\) and \(\beta\). We will also derive the score function as well.

1.1 Part A: Derivation of the Likelihood Function

To start, we need a sample. If \(X_1\), \(X_2\), \(...\), \(X_n\) are independent and identically distributed such that ~ \(\Gamma(\alpha, \beta)\), then we have the joint pdf: \[ f(x_i \mid \alpha, \beta) = \prod_{i=1}^{n} \frac{1}{\Gamma(\alpha)\beta^\alpha} x_i^{\alpha-1} e^{-x/\beta} \ \ \text{for} \ \ x > 0 \]

Simplifying this, we get our likelihood function: \[ L(\alpha, \beta) = [\frac{1}{\Gamma(\alpha)\beta^\alpha}]^n [\prod_{i=1}^{n}x_i]^{\alpha-1}[e^-[(1/\beta){\Sigma_{i=1}^nx_i}]] \]

Now, we take natural log of the likelihood function: \[\ln(L(\alpha, \beta) = n\]

2 Problem 2

For this problem, we are asked to estimate the parameters of a gamma distribution using a sample consisting of the lengths of time in hours spent by women in the delivery suite while giving birth (without a ceasarian section) for each of 7 days at John Radcliffe Hospital in Oxford, England. The data are taken from Davison (Statistical Models. Cambridge University Press, 2003).

birth <- c(2.1, 3.4, 4.25, 5.6, 6.4, 7.3, 8.5, 8.75, 8.9, 9.5, 9.75, 10, 10.4, 10.4, 16, 19, 4, 4.1, 5, 5.5, 5.7, 6.5, 7.25, 7.3, 7.5, 8.2, 8.5, 9.75, 11, 11.2, 15, 16.5, 2.6, 3.6, 3.6, 6.4, 6.8, 7.5, 7.5, 8.25, 8.5, 10.4, 10.75, 14.25, 14.5, 1.5, 4.7, 4.7, 7.2, 7.25, 8.1, 8.5, 9.2, 9.5, 10.7, 11.5, 2.5, 2.5, 3.4, 4.2, 5.9, 6.25, 7.3, 7.5, 7.8, 8.3, 8.3, 10.25, 12.9, 14.3, 4, 4, 5.25, 6.1, 6.5, 6.9, 7, 8.45, 9.25, 10.1, 0.2, 12.75, 14.6, 2, 2.7, 2.75, 3.4, 4.2, 4.3, 4.9, 6.25, 7, 9, 9.25, 10.7)

2.1 Part A: Maximum Likelihood Estimators

For this part, we will estimate the parameters via maximum likelihood method.

# gamma log-likelihood function
gamma.loglik <- function(params, data) {
  shape <- params[1]  # passing the parameters
  scale <- params[2]
  n <- length(birth)   # sample size
  
  # Log-likelihood for gamma distribution
loglik <- -n*log(gamma(shape)) - n*shape*log(scale) + 
  (shape - 1)*sum(log(data)) - sum(data)/scale
  
  return(loglik)
}

# Score (gradient) equation

gamma.score <- function(params, data) {
  shape <- params[1]
  scale <- params[2]
  n <- length(birth)
  
  # Gradient for shape parameter
  grad_shape <- -n*digamma(shape) - n*log(scale) + 
                sum(log(data)) 
  
  # Gradient for scale parameter
  grad_scale <- -(n*shape)/scale + sum(birth)/scale^2 
  
  return(c(grad_shape, grad_scale))
}


# Need to provide initial values for parameters
initial.params <- c(shape = 2.0, scale = 2.0)  # Reasonable starting values


# Using optim with Nelder-Mead method
mle.result <- optim(
  par = initial.params,
  fn = gamma.loglik,
  gr = gamma.score,
  data = birth,
  method = "L-BFGS-B",
  hessian = TRUE,
  control = list(trace = FALSE,
                 fnscale = -1,
                 maxit = 500,
                 abstol = 1e-8)
)
##
mle.result

$par
   shape    scale 
3.584610 2.125167 

$value
[1] -257.588

$counts
function gradient 
      14       14 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
          shape     scale
shape -30.53752 -44.70238
scale -44.70238 -75.40150

According to the output, our maximum likelihood estimators are \(\hat{\alpha} = 3.584613\) and \(\hat{\beta} = 2.125165\).

2.2 Part B: Method of Moments Estimators

This time, we will estimate the parameters via method of moments.

mom1 <- mean(birth)
mom2 <- mean(birth^2)

mm.est.scale <- mom1 / (mom2 - mom1^2)
mm.est.shape <- mom1 * mm.est.scale

mm.est.shape

[1] 4.424066

mm.est.scale

[1] 0.5807466

According to the output, our method of moments estimators are \(\tilde{\alpha} = 4.424066\) and \(\tilde{\beta} = 0.5807466\).

2.3 Part C: Which is Better?

So how do we compare these two methods of estimating parameters? For method of moments, we are following the idea that the expectation of a distribution is a function of the random variable. As for maximum likelihood estimation, we are following the idea that the entire distribution is a function of x. So we can work directly with the data/distribution when using maximum likelihood estimating instead of using means and variances or other moments. According to most sources on the internet, most seem to agree that maximum likelihood is much better and much more useful than method of moments. And based on the work i did for this assignment, i’m inclined to agree.