Maximum Likelihood Estimation (MLE)

1.) Gradient Derivation

For this analysis we will assume that our random variable X is Gamma distributed and the instances which make up our random sample of x resemble those distributive characteristics. The conditional PDF of X, whose values depend on the parameters of the Gamma Distribution - \(\alpha\) to represent shape and \(\beta\) to represent rate - is written below.

\[\Large f(x | \alpha, \beta) = \frac{1}{\Gamma(\alpha)\beta^\alpha}x^{(\alpha-1)}e^{(-x/\beta)} \hspace{0.5cm}x,\alpha, \beta> 0 \] To perform the partial derivations throughout this analysis, standard derivative rules such as the power, product and quotient rule were heavily relied upon. The quotient rule in particular and the manner in which I used it is depicted below.

\[Quotient \ Rule \ =\Large \frac{g(x)h'(x)-h(x)g'(x)}{[g(x)]^2}\]
\[\Large \frac{h(x) = x^{\alpha-1}e^{-x/\beta}}{g(x) = \Gamma(\alpha)\beta^\alpha} \]

Derivation on Original Function

To enhance my understanding of how \(\alpha\) and \(\beta\) individually impact the value of x in our distribution, I found the rate of change that each parameter is respectively responsible for (ie the partial derivative or score function).

Partial Derivative W/ Respect to Alpha

\[\Large 1: \frac{\partial f}{\partial \alpha} = e^{-x/\beta}* \frac{d}{d\alpha}[\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}] \]

\[\Large 2: \frac{d}{d\alpha}[\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}] = \frac{(\Gamma(\alpha)\beta^\alpha)[x^{\alpha-1}ln(x)]-(x^{\alpha-1})(\frac{d}{d\alpha}[\Gamma(\alpha)\beta^\alpha])}{[\Gamma(\alpha)\beta^\alpha)]^2}\]
\[\Large 3: \frac{d}{d\alpha}[\Gamma(\alpha)\beta^\alpha] = \Gamma(\alpha)*\beta^\alpha ln(\beta)+\beta^\alpha*\frac{d}{d\alpha}[\Gamma(\alpha)] \]
\[\Large 4: \frac{d}{d\alpha}[\Gamma(\alpha)\beta^\alpha] = \Gamma(\alpha)*\beta^\alpha ln(\beta)+\Gamma(\alpha)\psi^{(0)}(\alpha)*\beta^\alpha\]
\[\Large 5: \frac{d}{d\alpha}[\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}] = \frac{(\Gamma(\alpha)\beta^\alpha)[x^{\alpha-1}ln(x)]-(x^{\alpha-1})[\Gamma(\alpha)*\beta^\alpha ln(\beta)+\Gamma(\alpha)\psi^{(0)}(\alpha)*\beta^\alpha]}{[\Gamma(\alpha)\beta^\alpha)]^2}\]
\[\Large 6: \frac{d}{d\alpha}[\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}] = \frac{x^{\alpha-1}[ln(x)-ln(\beta)+\psi^{(0)}(\alpha)]}{\Gamma(\alpha)*\beta^\alpha} \]
Now multiply constant component with derivative component \[\Large 7: \frac{\partial f}{\partial \alpha} = e^{-x/\beta}* \frac{d}{d\alpha}[\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}] \]

\[\Large 8: \frac{\partial f}{\partial \alpha} = \frac{e^{-x/\beta}*x^{\alpha-1}[ln(x)-ln(\beta)+\psi^{(0)}(\alpha)]}{\Gamma(\alpha)*\beta^\alpha} \]
\[\Large 9: \frac{\partial f}{\partial \alpha} = \frac{x^{\alpha-1}e^{-x/\beta}}{\Gamma(\alpha)*\beta^\alpha} * [ln(x)-ln(\beta)+\psi^{(0)}(\alpha)] \]

\[\Large 10: \frac{\partial f}{\partial \alpha} = f(x|\alpha, \beta) * [ln(x)-ln(\beta)+\psi^{(0)}(\alpha)] \]

Partial Derivative W/ Respect to Beta

\[\Large 1: \frac{\partial f}{\partial \beta} = \frac{x^{\alpha-1}}{\Gamma(\alpha)} *\frac{d}{d\beta}[\frac{e^{(-x/\beta)}}{\beta^\alpha}] \]
Isolate derivative component \[\Large 2: \frac{d}{d\beta}[\frac{e^{-x/\beta}}\beta^\alpha] = \frac{\beta^\alpha * \frac{d}{d\beta}[e^{-x/\beta}]-e^{-x/\beta}*\frac{d}{d\beta}[\beta^\alpha]}{\beta^{2\alpha}} \]

\[\Large 3: \frac{d}{d\beta}[\frac{e^{-x/\beta}}\beta^\alpha] = \frac{\beta^\alpha*\frac{-x}{\beta^2}*e^{-x/\beta}-e^{-x/\beta}*\alpha\beta^{\alpha-1}}{\beta^{2\alpha}} \]

\[\Large 4: \frac{d}{d\beta}[\frac{e^{-x/\beta}}{\beta^\alpha}] = \frac{(\beta^{\alpha-2} * x*e^{-x/\beta})-(e^{-x/\beta}*\alpha\beta^{\alpha-1})}{\beta^{2\alpha}} \]

Now multiply constant component with derivative component \[\Large 5: \frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})[\beta^{\alpha-2}xe^{-x/\beta}-e^{-x/\beta}\alpha\beta^{\alpha-1}]}{\Gamma(\alpha)\beta^{2\alpha}}\]
\[\Large 7: \frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})(e^{-x/\beta})(\beta^{\alpha-1})(\beta^{-1}x-\alpha)}{\Gamma(\alpha)\beta^{2\alpha}} \]
\[\Large 8: \frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})(e^{-x/\beta})(\beta^{-1}x-\alpha)}{\Gamma(\alpha)\beta^{\alpha+1}} \]
\[\Large 9:\frac{\partial f}{\partial \beta} =\frac{(x^{\alpha-1})(e^{-x/\beta})(\beta^{-1}x-\alpha)(\beta^{-(\alpha+1)})}{\Gamma(\alpha)} \]
\[\Large 10:\frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})(e^{-x/\beta})}{\Gamma(\alpha)}*\frac{(\beta^{-1}x-\alpha)*\beta^{-\alpha}*\beta^{-1}}{1} \]
\[\Large 11: \frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})(e^{-x/\beta})}{\Gamma(\alpha)\beta^\alpha}*\frac{(\beta^{-1}x-\alpha)*\beta^{-1}}{1} \]
\[\Large 12: \frac{\partial f}{\partial \beta} = \frac{(x^{\alpha-1})(e^{-x/\beta})}{\Gamma(\alpha)\beta^\alpha}*\frac{\beta^{-2}x-\beta^{-1}\alpha}{1} \]
\[\Large 13: \frac{\partial f}{\partial \beta} = {\frac{(x^{\alpha-1})(e^{-x/\beta})}{\Gamma(\alpha)\beta^\alpha}}*[\frac{x}{\beta^2}-\frac{\alpha}{\beta}]\]
\[\Large 14: \frac{\partial f}{\partial \beta} =f(x | \alpha, \beta) *[\frac{x}{\beta^2}-\frac{\alpha}{\beta}] \]

Gradient Vector of Original Function

Each of these partial derivatives then compose the gradient vector of the original density function. The top component \(\frac{\partial f}{\partial \alpha}\) tell us how a change in alpha’s value will affect the value of x, assuming all other parameters (like \(\beta\)) remain constant. Additionally, the bottom component tell us analogous information regarding \(\beta\).

\[\Large \nabla f_{\alpha,\beta} = \begin{bmatrix} \frac{\partial f}{\partial \alpha} \\ \frac{\partial f}{\partial \beta} \end{bmatrix} = \begin{bmatrix} f(x|\alpha, \beta) * [ln(x)-ln(\beta)+\psi^{(0)}(\alpha)] \\ f(x | \alpha, \beta) *[\frac{x}{\beta^2}-\frac{\alpha}{\beta}] \end{bmatrix} \]

1a.) Log Likelihood Transformation

The likelihood function tells us how probable or ‘likely’ our parameter of interest (usually denoted as \(\theta\)) is equal to what we believe it is equal to given our observed data. It is mathematically defined as the product of each observation’s conditional probability operating under the assumption that \(\theta\) = \(\theta\). In this instance, since the value of x is conditional upon the values of \(\alpha\) and \(\beta\), the likelihood function is the following:

\[\Large L(\theta | x) = \prod_{i=1}^{n} f(x_i | \theta) = \prod_{i=1}^{n} f(\alpha_i, \beta_i | \theta). \]
Given that product series can be notoriously difficult to work with, the log-likelihood function is more often used for parameter estimation as logs transform multiplicative procedures into additional ones. The defining of the log-likelihood function for our original given PDF is below.

\[\Large 1: ln[L(\theta|x)] = \sum_{i=1}^{n} ln[f(x_i | \theta)] = \sum_{i=1}^{n} ln[f(\alpha_i, \beta_i | \theta)] \]

\[\Large 2: ln[f(\alpha_i, \beta_i | \theta)] = ln[x^{\alpha-1}e^{-x/\beta}]-ln[\Gamma(\alpha)\beta^\alpha]\]

\[\Large 3: ln[f(\alpha_i, \beta_i | \theta)] = ln[x^{\alpha-1}]+ln[e^{-x/\beta}]-ln[\Gamma(\alpha)]-ln[\beta^\alpha]\]

\[\Large 4: ln[f(\alpha_i, \beta_i | \theta)] = (\alpha-1)ln(x) - \frac{x}{\beta}-ln[\Gamma(\alpha)]-\alpha ln[\beta]\]

\[\Large 5: ln[L(\theta | x)] = \sum_{i=1}^{n} (\alpha_i-1)ln(x_i) - \frac{x_i}{\beta_i}-ln[\Gamma(\alpha_i)]-\alpha_i ln[\beta_i]\]

1b.) Log Likelihood Gradient Derivation

Just as I previously found the partial derivatives of the density function, now I will compose the gradient for the log-likelihood function. This will tell us how a change in either parameter of interest - \(\alpha\) or \(\beta\) - affects the log-likelihood of our targeted estimation parameter \(\theta\). Calculations and breakdowns are below.

Partial LOG Derivative W/ Respect to Alpha

\[\Large 1: ln[\frac{\partial f}{\partial \alpha}] = \frac{d}{d\alpha}[\alpha ln(x) - ln(x)] - \frac{d}{d\alpha}[ln[\Gamma(\alpha)]] - ln[\beta] \]

\[\Large 2: ln[\frac{\partial f}{\partial \alpha}] = ln(x) - ln[\beta] - \frac{d}{d\alpha}[ln[\Gamma(\alpha)]]\]

\[\Large 3: ln[\frac{\partial f}{\partial \alpha}] = ln(x) - ln[\beta] - \psi^{(0)}(\alpha)]\]

Partial LOG Derivative W/ Respect to Beta

\[\Large 1: ln[\frac{\partial f}{\partial \beta}] = \frac{d}{d\beta}[\frac{-x}{\beta} - \alpha ln[\beta]]\]

\[\Large 2: ln[\frac{\partial f}{\partial \beta}] = x\beta^{-1-1}-\alpha*\frac{1}{\beta}\]

\[\Large 3: ln[\frac{\partial f}{\partial \beta}] = \frac{x}{\beta^2}-\frac{\alpha}{\beta}\]

Gradient Vector of LOG Likelihood Function

\[\Large \nabla ln[f_{\alpha,\beta}] = \begin{bmatrix} ln[\frac{\partial f}{\partial \alpha}] \\ ln[\frac{\partial f}{\partial \beta}] \end{bmatrix} = \begin{bmatrix} [ln(x_i)-ln(\beta)+\psi^{(0)}(\alpha)] \\ [\frac{x_i}{\beta^2}-\frac{\alpha}{\beta}] \end{bmatrix} \]

2.) Birth Data

The available data we are examining is a sample of 95 womens’ time spent (hours) in the delivery suite while giving birth at John Radcliffe Hospital in Oxford, England. The study was conducted decades ago. For greater context, summary statistics and a histogram visualization of the sample distribution are below.

Births = 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)

Births_Data = data.frame(Births)
sum_stats_table = function(variable){
  table = cbind(
    round(min(variable),2),
    round(quantile(variable, 0.25),2),
    round(median(variable),2),
    round(mean(variable),2),
    round(quantile(variable, 0.75),2),
    round(max(variable),2)
  )
  colnames(table) = c("Min", "Q1", "Med", "Mean", "Q3", "Max")
  rownames(table) = NULL
  return(table)
}
  # I created a table for summary statistics of a variable, can use this moving forward instead of manually putting the summary stats in table form


Summary_Table = sum_stats_table(variable = Births)

kable(Summary_Table,align = "c",
      caption = "<span style='color:##000000;'>Summary Stats of Birth Times (Hours)</span>") %>%
  kable_styling(
    bootstrap_options = c("striped", "bordered"),
    full_width = FALSE,
    position = "center")

Summary Stats of Birth Times (Hours)

Min	Q1	Med	Mean	Q3	Max
1.5	4.95	7.5	7.72	9.75	19

hist(Births, 
     main = "Distribution of Delivery Times",
     xlab = "Delivery Time (Hours)",
     ylab = "Count",
)
box(lwd=3)

2a.) Using MLE to Approximate Alpha and Beta

Operating under the assumption that the broader population of delivery times for women at said hospital at said time period is gamma distributed, we will utilize the log-likelihood gradient derived earlier to estimate \(\alpha\) and \(\beta\). These estimations will be denoted as \(\hat{\alpha}\) and \(\hat{\beta}\). Note that the process of MLE is simply setting a partial derivative equal to zero, and solving for the isolated parameter of interest (first \(\alpha\), and then \(\beta\) in this case).

## Optim Estimation Redo

 # note: in R, there is no ln function, but the log function defaults to natural log

# Code is modeled after section 5.2.2 in class notes

# Overall log likelihood
Gamma_LogLikelihood = function(parameters, data){
  alpha = parameters[1]
  beta = parameters[2]
  
  n = length(data)
  
  LogLike = (n * alpha * log(beta)) - 
    (n * log(gamma(alpha))) +
    (alpha - 1) * sum(log(data)) - 
       (beta * sum(data))
}
    # I originally made programming error of having x as a parameter in my function instead of simply summing the data


# Gradient Equations
Gamma_Gradient = function(parameters, data){
  alpha = parameters[1]
  beta = parameters[2]
  
  n = length(data)
  
  # Alpha Gradient
  alpha_grad = sum(log(data)) - (n*log(beta)) + (n*digamma(alpha))
    # Simply distribute the n through and use the sum function for all iterations of the "x" datapoint in our birthing dataset
  
  # Beta Gradient
  beta_grad = ((sum(data)/(beta^2)) - (n*(alpha/beta)))
  
  return(c(alpha_grad, beta_grad))  
}

Starting_Parameters = c(alpha = 3, beta = 3)
  # Arbitrarily setting starting parameters

MLE = optim(
  par = Starting_Parameters,
  fn = Gamma_LogLikelihood,
  gr = Gamma_Gradient,
  data = Births,
#  method = "L-BFGS-B",
  hessian = TRUE,
  control = list(trace = FALSE,
                 fnscale = -1,
                 maxit = 500,
                 abstol = 1e-8)
)

invisible(
  round(MLE$par[1],3) # alpha 
)

invisible(
  round(MLE$par[2],3) # beta 
)

Via R’s optim() function, I found an approximate value of 4.389 for \(\hat{\alpha}\) (our shape parameter) and 0.568 for \(\hat{\beta}\) (our rate parameter). For further confirmation of reasonable accuracy, I plotted a Gamma distribution with said parameter values alongside the distribution of the birthing times from our sample data.

hist(Births, 
     main = "Distribution of Delivery Times",
     xlab = "Delivery Time (Hours)",
     ylab = "Probability Density",
     ylim = c(0, 0.13),
     prob = TRUE
)

# alpha = shape
# beta = scale
curve(dgamma(x , shape = 4.3885430, rate = 0.5681101), 
      col = "red",
      lwd = 3,
      add = TRUE)

legend("topright",
       legend = "Gamma (4.389, 0.568)",
       col = "red",
       lwd = 3,
       box.lwd=3)

box(lwd=3)

2b.) Using Moment Estimation to Approximate Alpha and Beta

Alongside using MLE, I chose to additionally approximate \(\alpha\) and \(\beta\) via moment estimation. The estimated parameters found via this method will be denoted as \({\widetilde\alpha}\) and \({\widetilde\beta}\).

## kth moment of distribution = summing datasets values raised to kth power, and then dividing by n

moment1 = mean(Births) # = 7.723158
moment2 = mean(Births^2) # = 72.24447
variance = moment2 - (moment1^2) # = 12.59731

# mean of Gamma distribution is alpha/beta
# var of Gamma distribution is alpha/(beta^2)
  ### Formulas are for shape parameter alpha and rate parameter beta, there are other versions of Gamma notation

Per the general formula of moment estimation, the kth moment occurs when all of a datasets’ values are raised to the kth power, summed and then divided by the sample size, n. Combining that with the expected value and variance formulas specific to the Gamma distribution, we can create the system of equations below.

\[\Large Moment 1 = \frac{\alpha}{\beta} = \bar{x} \approx 7.723 \] \[\Large Moment 2 = \frac{\alpha^2}{\beta^2} = \bar{x^2} \approx 72.244 \]

\[\Large Var = \frac{\alpha}{\beta^2} = \bar{x^2} - [\bar{x}]^2 \approx 12.597 \]

Now using this system of equations, algebraic manipulation is performed to isolate and estimate the values of our parameters \(\alpha\) and \(\beta\).

\[\Large \frac{\alpha}{\beta^2} = \frac{1}{\beta} * \frac{\alpha}{\beta}\]

\[\Large \frac{var}{mean} = \frac{1}{\beta} \]

\[\Large \beta = \frac{\bar{x}}{\bar{x^2} - [\bar{x}]^2} \]

\[\Large \alpha = \bar{x} * \beta \]

mom_beta = moment1/variance
mom_alpha = moment1*mom_beta

invisible(
  round(mom_alpha,3)
)

invisible(
  round(mom_beta,3)
)

After performing all calculations, the moment estimations are the following; \({\widetilde\alpha}\) = 4.735 and \({\widetilde\beta}\) = 0.613. These estimates are very similar to those found via MLE in the previous section. For enhanced comparison, gamma curves with adherence to each estimation method’s outcomes are overlaid on top of our sample data below.

hist(Births, 
     main = "Distribution of Delivery Times",
     xlab = "Delivery Time (Hours)",
     ylab = "Probability Density",
     ylim = c(0, 0.13),
     prob = TRUE
)

# alpha = shape
# beta = scale
curve(dgamma(x , shape = 4.3885430, rate = 0.5681101), 
      col = "red",
      lwd = 3,
      add = TRUE)

curve(dgamma(x , shape = 4.734915, rate = 0.6130801), 
      col = "blue",
      lwd = 3,
      add = TRUE)

legend("topright",
       legend = c("MLE (4.389, 0.568)", "MOM (4.735, 0.613)"),
       col = c("red", "blue"),
       lwd = 3,
       box.lwd=3)

box(lwd=3)

2c.) Moment Estimation vs MLE

Broadly speaking, both moment estimation (MME) and maximum likelihood estimation (MLE) are fundamentally sound methods of estimating population parameters. MME consists of creating a system of equations with heavy reliance on the kth moment principle. MLE on the other hand involves finding partial derivatives of the log likelihood function, equating them to zero and calculating which value of \(\theta\) makes our sample data most probable.

With regards to accuracy, MLE is generally the preferable approach, yielding results with lower variance and less risk of biased influence. Moment estimation can be particularly vulnerable to bias when sample sizes are not sufficiently large. On the other hand, likelihood estimation requires we know the applicable density function and that the partial derivative of that function can be reasonably found.

To summarize, resorting to moment estimation is ideal when we do now know a population’s density function or said function’s partial derivative goes beyond standard calculus-based/numerical techniques. However, if the density function is readily available, then likelihood estimation is preferred due to its greater accuracy.