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.

I will begin by finding the likelihood function.

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

This likelihood function can be rewritten as follows.

\(L(\alpha, \beta) = \left[\frac{1}{\Gamma(\alpha)\beta^\alpha}\right]^n \prod_{i=1}^{n} x_i^{\alpha-1} \exp\left(-\frac{1}{\beta}\sum_{i=1}^{n} x_i\right).\)

Next, I will take the logarithm of this likelihood function in order to obtain the log-likelihood function.

\(\ell(\alpha, \beta) = \log L(\alpha, \beta)\)

Taking the logarithm of the likelihood function results in:

\(\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\)

Therefore, this is the log-likelihood function of the Gamma distribution with parameters \(\alpha\) and \(\beta\).

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

First, I will derive the score function for \(\alpha\) from the log-likelihood function that I found in the previous part. I will take the partial derivative with respect to the \(\alpha\) parameter. This will give the score function.

To begin, I know that \(\frac{d}{d\alpha}\log \Gamma(\alpha) = \psi_0(\alpha)\) . So, I can use this fact to help with taking the partial derivative of the log-likelihood function with respect to the \(\alpha\) parameter. Doing so results in the following:

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

Therefore the score function for \(\alpha\) is: \(s_\alpha(\alpha,\beta) = \sum_{i=1}^{n} \log x_i - n\log \beta - n\psi(\alpha)\)

Now, I will find the score function for \(\beta\). I will begin by taking the derivative with respect to \(\beta\) to get the partial derivative for the \(\beta\) parameter. This will give the score function.

To begin, I will differentiate the log-likelihood function with respect to the \(\beta\) parameter. From the original log-likelihood functions, there are two terms with the \(\beta\) parameter. I know these derivatives are as follows, \(\frac{\partial}{\partial \beta} (-n\alpha \log \beta) = -\frac{n\alpha}{\beta}\) and \(\frac{\partial}{\partial \beta} \left(-\frac{1}{\beta}\sum_{i=1}^{n} x_i\right) = \frac{1}{\beta^2}\sum_{i=1}^{n} x_i\)

So, the partial derivative of the log-likelihood function altogether results in the following.

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

Therefore, the score function for \(\beta\) is: \(s_\beta(\alpha,\beta) = -\frac{n\alpha}{\beta} + \frac{1}{\beta^2}\sum_{i=1}^{n} x_i\)

So altogether, the score functions for \(\alpha\) and \(\beta\), respectively are the following:

\(s_\alpha(\alpha,\beta) = \sum_{i=1}^{n} \log x_i - n\log \beta - n\psi(\alpha)\)

\(s_\beta(\alpha,\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.

I will begin by storing the data values in a data set that will be named ‘time’.

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

Next, I will find the MLEs for \(\alpha\) and \(\beta\) by using the optim() function and the score functions I found in question #1 part b.

n <- length(time)

# Log-likelihood function
loglik <- function(par){
  alpha <- par[1]
  beta  <- par[2]
  
  if(alpha <= 0 || beta <= 0) return(Inf)
  
  -(-n*lgamma(alpha) 
    - n*alpha*log(beta)
    + (alpha-1)*sum(log(time))
    - sum(time)/beta)}

# Gradient functions from question #1
grad <- function(par){
  alpha <- par[1]
  beta  <- par[2]
  
  d_alpha <- sum(log(time)) - n*log(beta) - n*digamma(alpha)
  d_beta  <- -n*alpha/beta + sum(time)/beta^2
  
  return(-c(d_alpha, d_beta))
}

mean_time <- mean(time)
var_time <- var(time)

alpha_init <- mean_time^2 / var_time
beta_init  <- var_time / mean_time

# Using the optim() function to get the MLEs
optim(c(alpha_init, beta_init),
      loglik,
      grad,
      method="L-BFGS-B",
      lower=c(1e-6,1e-6))

$par
[1] 4.387855 1.760122

$value
[1] 251.1173

$counts
function gradient 
       8        8 

$convergence
[1] 0

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

It turns out that \(\hat{\alpha}\) = 4.388 and \(\hat{\beta}\) = 1.760.

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

In this part, I will apply the method of moments to obtain the estimators for \(\alpha\) and \(\beta\). For a Gamma distribution, I know that \(E(X) = \alpha \beta\) and Var(X) = \(\alpha \beta^2\). I will write these in terms of the sample mean \(\bar{x}\) and the sample variance \(s^2 = \alpha \beta^2\).

\(\tilde{\alpha} = \frac{\bar{x}^2}{s^2}\)

\(\tilde{\beta} = \frac{s^2}{\bar{x}}\)

Now, to calculate the values of \(\tilde{\alpha}\) and \(\tilde{\beta}\) based upon these moments.

# Method of moments
alpha_mom <- mean_time^2 / var_time
beta_mom  <- var_time / mean_time

alpha_mom

[1] 4.685073

beta_mom

[1] 1.64846

It turns out that \(\tilde{\alpha}\) = 4.685 and \(\tilde{\beta}\) = 1.648.

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.

I used both the method of moments estimation and the maximum likelihood estimation to estimate the parameters \(\alpha\) and \(\beta\). Both the MLE and the MME methods resulted in similar values. The MLE method resulted in \(\hat{\alpha}\) = 4.388 and \(\hat{\beta}\) = 1.760. The MME method resulted in \(\tilde{\alpha}\) = 4.685 and \(\tilde{\beta}\) = 1.648. These values are quite similar showing that both methods provide use in estimating the parameters.

A major advantage of the Method of moments estimation method was that it is simpler and quicker to compute. For this method, I used facts I knew about the Gamma distribution and its mean and variance in order to get the sample moments, and with some short calculations, this gave me my parameter estimates. On the other hand, the maximum likelihood function required some more complex calculations which required the use of the R optim() function. This required a few extra steps and work to obtain the MLE estimations due to the additional computations required, and needing the use of the gradient functions I calculated in question #1 for the partial derivatives to get the score functions.

However, the MLE method does bring some advantages of its own. The MLE has strong asymptotic normality which brings consistency and states that the distribution of the estimators, in this case \(\alpha\) and \(\beta\) will approach a Normal distribution as the sample size increases to infinity. This property makes the MLE method consistent, reliable, and efficient.

Some limitations of the MME method include that it can be prone to bias, and this can happen especially in skewed distributions. This could occur if the sample moments are not the best representations of the true parameters they are estimating. Therefore, the MME does have some drawbacks and limitations when it comes to potential bias in its estimations if the population is skewed and therefore the sample moments do not provide the greatest accuracy towards the true parameters. On the other hand, the MLE method can also have some drawbacks. While the MLE method is great for large samples, in the case of small samples it can have some potential for bias as well. Additionally, the MLE method requires more complexity in its calculations especially when compared to the MME method.

Overall, I would use the Maximum likelihood estimation as the preferred approach due to its consistency and reliability, even for large sample sizes. The sample size in this scenario was significantly large, meaning that the issue of the MLE having potential for bias in small sample sizes will not be relevant in this case due to a large sample size. Additionally, even though the MLE method requires more complex calculations than the MME method, this trade off is worthwhile due to the consistency and reliability of its calculations. The MLE method is efficient and reliable, making it a great choice for estimation. In this case, I would recommend the MLE method for these reasons.