Assignment 4: 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.

Policies of Using AI Tools

Policy on AI Tool Use: Students must adhere to the AI tool policy specified in the course syllabus. The direct copying of AI-generated content is strictly prohibited. All submitted work must reflect your own understanding; where external tools are consulted, content must be thoroughly rephrased and synthesized in your own words.

Code Inclusion Requirement: Any code included in your essay must be properly commented to explain the purpose and/or expected output of key code lines. Submitting AI-generated code without meaningful, student-added comments will not be accepted.

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, \]

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 log-likelihood function with respect to the gamma distribution parameters \(\alpha\) and \(\beta\). To this end,

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

Answer to Question 1 Part A

The Gamma distribution is defined by the probability density function:

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

where \(\alpha > 0\) is the shape parameter and \(\beta>0\) is the scale parameter.

Therefore the likelihood function is:

\[ L(\alpha,\beta \mid x)=\prod_{i=1}^{n}{f(x_i\mid\alpha,\beta)}=\prod_{i=1}^{n}{\frac{1}{\Gamma(\alpha)\beta^{\alpha}}x_i^{\alpha-1}e^{-x_i/\beta}} \]

And the log-likelihood function is:

\[ \ell(\alpha,\beta\mid x)=\ln{L(\alpha,\beta \mid x)}=\sum_{i=1}^{n}{\ln{\left[\frac{1}{\Gamma(\alpha)\beta^{\alpha}}x_i^{\alpha-1}e^{-x_i/\beta}\right]}} \]

This can be further simplified to:

$$ \[\begin{aligned} \ell(\alpha,\beta\mid x) &= \sum_{i=1}^{n}{\ln{\left[\frac{1}{\Gamma(\alpha)\beta^{\alpha}}x_i^{\alpha-1}e^{-x_i/\beta}\right]}} \\\\ &= \sum_{i=1}^{n}{\left[\ln{1}-(\ln{\Gamma(\alpha)}+\alpha\ln{\beta})+(\alpha-1)\ln{x_i}-\frac{x_i}{\beta}\right]}\\\\ &= -n\ln{\Gamma(\alpha)}-n\alpha\ln{\beta}+(\alpha-1)\sum_{i=1}^{n}{\ln{x_i}}-\frac{1}{\beta}\sum_{i=1}^{n}{x_i} \end{aligned}\]

$$

Therefore, the simplifed log-liklihood function is:

\[ \ell(\alpha,\beta\mid x)=-n\ln{\Gamma(\alpha)}-n\alpha\ln{\beta}+(\alpha-1)\sum_{i=1}^{n}{\ln{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 log-likelihood function in part (a).

Answer to Question 1 Part B

To find the first score function we will differentiate with respect \(\beta\).

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\beta} &= \frac{\partial}{\partial\beta}\left[-n\ln{\Gamma(\alpha)}-n\alpha\ln{\beta}+(\alpha-1)\sum_{i=1}^{n}{\ln{x_i}}-\frac{1}{\beta}\sum_{i=1}^{n}{x_i}\right]\\\\ &= -n\alpha\frac{1}{\beta}+\frac{1}{\beta^2}\sum_{i=1}^{n}{x_i} \end{aligned}\]

$$

To find the second score function we will differentiate with respect of \(\alpha\):

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\alpha} &= \frac{\partial}{\partial\alpha} \left[-n\ln{\Gamma(\alpha)}-n\alpha\ln{\beta}+(\alpha-1)\sum_{i=1}^{n}{\ln{x_i}}-\frac{1}{\beta}\sum_{i=1}^{n}{x_i}\right] \\\\ &= -n\frac{\Gamma (\alpha)\psi_0(\alpha)}{\Gamma(\alpha)}-n\ln{\beta}+\sum_{i=1}^{n}{\ln{x_i}} \end{aligned}\]

$$

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.

Answer to Question 2 Part A

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, 
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)

#Used to enter log-logistic function
gamma.loglik <- function(params, data) {
  shape <- params[1]  # Used to get initial parameters
  scale <- params[2]
  n <- length(data)   # Used to get sample size
  

  loglik <- -n*log(gamma(shape))-n*shape*log(scale)+(shape-1)*sum(log(data))-(1/scale)*sum(data) 
  
  return(loglik)    # Return negative for minimization
}


# Used to enter gradient functions

gamma.score <- function(params, data) {
  shape <- params[1] #Used to get initial parameters
  scale <- params[2]
  n <- length(data) #Used to get sample size
  
  # Gradient for shape parameter
  grad_shape <-  -n*digamma(shape)-n*log(scale)+sum(log(data))
  
  # Gradient for scale parameter
  grad_scale <- -n*shape*(1/scale)+(1/(scale**2))*sum(data)
  
  return(c(grad_shape, grad_scale))
}


# Entering initial guessed parameters
initial.params <- c(shape = 2.2, scale = 1.5)  


# Using optim and previous functions to get results
mle.result <- optim(
  par = initial.params,
  fn = gamma.loglik,
  gr = gamma.score,
  data = birth,
  method = "L-BFGS-B",
  lower = c(1e-6, 1e-6),
  hessian = TRUE,
  control = list(trace = FALSE,
                 fnscale = -1,
                 maxit = 500,
                 abstol = 1e-8)
)
##
mle.result #printing mle.result determined earlier

$par
   shape    scale 
4.387881 1.760111 

$value
[1] -251.1173

$counts
function gradient 
      15       15 

$convergence
[1] 0

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

$hessian
          shape      scale
shape -24.30316  -53.97386
scale -53.97386 -134.55458

Therefore the MLE estimations for \(\alpha\) and \(\beta\) are \(\hat{\alpha}=4.387881\) and \(\hat{\beta}=1.760111\).

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

Answer to Question 2 Part B

As determined in the 5th Week Lecture, the method of moment estimators for a gamma distribution are:

$$ \[\begin{aligned} \hat{\mu} &= \frac{1}{n}\sum_{i=1}^{n}{X_i}=\bar{X}\\\\ \hat{\sigma^2} &= \frac{1}{n}\sum_{i=1}^{n}{X_i^2-\bar{X}^2}=\frac{1}{n}\sum_{i=1}^{n}{(X_i-\bar{X})^2} \end{aligned}\]

$$

We can use these method of moment estimators to obtain estimators for \(\alpha\) and \(\beta\) for our datset:

# Method of Moments estimation
mom_mu <- mean(birth) #Used to get moments from sample
mom_sigma_sq <- mean(birth^2) - mom_mu^2
mom_sigma <- sqrt(mom_sigma_sq)

cat("MoM estimates: μ_tilde =", round(mom_mu, 3), #Used to print moments
    ", σ_tilde =", round(mom_sigma, 3), "\n")

MoM estimates: μ_tilde = 7.723 , σ_tilde = 3.549 

Therefore the MME estimations for \(\alpha\) and \(\beta\) are \(\tilde{\alpha}=7.723\) and \(\tilde{\beta}=3.549\).

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 to Question 2 Part C

Method of moments estimation and maximum likelihood estimation are both methods of estimating the parameters of a distribution using sample data. Method of moments attempts to estimate these parameters by matching moments of the sample to moments of a theoretical probability model. Alternatively, maximum likelihood estimation considers how likely certain parameter values are given the sampled data. Given that MME only considers a few specific moments, it is less powerful, but also less computationally intense. Alternatively, MLE is more computationally intense but more powerful method of estimating parameters. MME would most likely be better if you just need a rough estimate of the parameters, while MLE might be better if you need a result that takes into account the entire shape of a probability distribution.