Assignment 5: Maximum Likelihood Estimation

Assignment Objectives

• Comprehend the likelihood function and its properties.

• Master the maximum likelihood estimation framework and required components.

• Understand the plug-in principle underlying MLE.

• Implement maximum likelihood estimation procedures in R.

Gamma Distribution Revisited

Let \(X\) be the two parameter Gamma random variable with density function

\[ f(x \mid \alpha, \beta) = \frac{1}{\Gamma(\alpha)\beta^\alpha} x^{\alpha-1} e^{-x/\beta} \ \ \text{for} \ \ x > 0, \] L(α,β∣x1​,…,xn​)=i=1∏n​Γ(α)βα1​xiα−1​e−xi​/β.

where with \(\alpha > 0\) (shape), \(\beta>0\) (scale), and

\[ \Gamma(\alpha) = \int_{0}^{\infty} t^{\alpha-1} e^{-t} \, dt, \quad \alpha > 0. \]

\(\Gamma(x)\) can be computed in R using the gamma() function. The derivative of \(\Gamma(x)\) are given respectively by

\[ \Gamma^\prime(z) = \Gamma (z)\psi_0(z) \]

where \(\psi_0(z) = \frac{d}{dz} \ln \Gamma{z}\). In R, the digamma function \(\psi_0(z)\) is evaluated by digamma().

This assignment focuses on finding maximum likelihood estimate of parameters \(\alpha\) and \(\beta\) based on a real-world application data set.

Question 1: Derive gradient (first order partial derivative) of likelihood function

Assume that \(\{x_1, x_2, \cdots, x_n \} \to \text{ gamma }(\alpha, \beta)\) with density function given by

\[ f(x \mid \alpha, \beta) = \frac{1}{\Gamma(\alpha)\beta^\alpha} x^{\alpha-1} e^{-x/\beta} \ \ \text{for} \ \ x > 0, \]

Derive the gradient of the likelihood function with respect to the gamma distribution parameters \(\alpha\) and \(\beta\). To this end,

a). Write out the full likelihood function based on the given data and the density function provided above.

The likelihood function is \[ L(\alpha,\beta \mid x_1,\ldots,x_n)=\prod_{i=1}^{n}\frac{1}{\Gamma(\alpha)\beta^{\alpha}}\, x_i^{\alpha-1}\, e^{-x_i/\beta}. \] Taking logs, the log-likelihood is \[ \ell(\alpha,\beta)=- n\log\Gamma(\alpha)- n\alpha\log\beta+(\alpha-1)\sum_{i=1}^{n}\log x_i- \frac{1}{\beta}\sum_{i=1}^{n} x_i. \]

b). Derive the score functions (the gradient of the log-likelihood) from the full likelihood function in part (a).

Using \(\dfrac{d}{d\alpha}\log\Gamma(\alpha)=\psi_0(\alpha)\), the score for \(\alpha\) is

\[ \frac{\partial \ell(\alpha,\beta)}{\partial \alpha} = -n\psi_0(\alpha) -n\log\beta +\sum_{i=1}^{n}\log x_i. \]

The score for \(\beta\) is

\[ \frac{\partial \ell(\alpha,\beta)}{\partial \beta} = -\frac{n\alpha}{\beta} +\frac{1}{\beta^2}\sum_{i=1}^{n}x_i. \]

Question 2: Birth data set

The following R code reads in a data set containing, for each of 7 days, the lengths of time in hours spent by women in the delivery suite while giving birth (without a ceasarian section) at John Radcliffe Hospital in Oxford, England. The data are taken from Davison (Statistical Models. Cambridge University Press, 2003).

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, 
10.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

Assume the data are generated from a gamma distribution. The objective is to use these data and the designated algorithm to find the maximum likelihood estimates (MLEs) of the parameters \(\alpha\) and \(\beta\).

a). Find the MLEs of \(\alpha\) and \(\beta\), denoted by \(\hat{\alpha}\) and \(\hat{\beta}\), using gradient-based optimization via the R function optim() with the gradient vector derived in Question 1.

# Birth time data (hours)
x <- 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,
  10.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
)

n <- length(x)
sx <- sum(x)
slogx <- sum(log(x))
xbar <- mean(x)
s2 <- var(x)   # sample variance (uses n-1 denominator)

# Negative log-likelihood for Gamma(shape=alpha, scale=beta)
negloglik_gamma <- function(par, x) {
  alpha <- par[1]
  beta  <- par[2]
  if (alpha <= 0 || beta <= 0) return(Inf)

  n <- length(x)
  ll <- -n * lgamma(alpha) - n * alpha * log(beta) + (alpha - 1) * sum(log(x)) - (1 / beta) * sum(x)
  return(-ll)  # negative for minimization
}

# Gradient of negative log-likelihood ( = - score )
neggrad_gamma <- function(par, x) {
  alpha <- par[1]
  beta  <- par[2]
  n <- length(x)

  dlda <- -n * digamma(alpha) - n * log(beta) + sum(log(x))
  dldb <- -(n * alpha) / beta + sum(x) / (beta^2)

  return(c(-dlda, -dldb))
}

# Starting values: use Method of Moments (good, stable starting point)
alpha_start <- xbar^2 / s2
beta_start  <- s2 / xbar
start <- c(alpha_start, beta_start)

mle_fit <- optim(
  par = start,
  fn  = negloglik_gamma,
  gr  = neggrad_gamma,
  x   = x,
  method = "L-BFGS-B",
  lower = c(1e-8, 1e-8)   # enforce positivity
)

alpha_hat <- mle_fit$par[1]
beta_hat  <- mle_fit$par[2]

data.frame(
  alpha_hat = alpha_hat,
  beta_hat  = beta_hat,
  convergence = mle_fit$convergence
)

  alpha_hat beta_hat convergence
1  4.387855 1.760122           0

b). Apply the method of moments to obtain estimators for \(\alpha\) and \(\beta\). Denote these moment estimators as \(\tilde{\alpha}\) and \(\tilde{\beta}\).

alpha_tilde <- xbar^2 / s2
beta_tilde  <- s2 / xbar

data.frame(alpha_tilde = alpha_tilde, beta_tilde = beta_tilde)

  alpha_tilde beta_tilde
1    4.685073    1.64846

c). Conduct a brief literature review comparing the method of moments estimation (MME) and maximum likelihood estimation (MLE). Synthesize the key advantages and limitations of each, concluding with a practical recommendation.

Answer: The Method of Moments (MME) estimates parameters by matching sample moments (like the mean and variance) to the theoretical moments of the distribution. It is simple to compute and easy to understand, but it may not always provide the most efficient estimates.

Maximum Likelihood Estimation (MLE) chooses the parameter values that make the observed data most likely under the assumed model. When the model is appropriate and the sample size is large, MLE typically produces more accurate and statistically efficient estimates. However, it often requires numerical optimization.

In practice, MME is useful for quick estimation and as starting values, while MLE is generally preferred for final inference due to its stronger theoretical properties.