Assignment 3: Methods of Moment Estimation

Assignment Objectives

• Master the fundamental concepts of point estimation and performance metrics

• Understand the theoretical foundation of the method of moments estimator (MME)

• Implement MME in R, incorporating numerical approximation methods

Use of 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.

Log-logistic Distribution

The log-logistic distribution (also known as the Fisk distribution) is a continuous probability distribution that is particularly useful in contexts where data exhibit non-negative, skewed behavior and where the hazard rate is unimodal (increases to a peak and then decreases). It has been widely used in the areas such as survival analysis and reliability engineering, environmental science, economics, pharmacology, finance and risk management, etc.

For given shape parameter \(\beta\) and scale parameter \(\alpha\), the cumulative distribution function

\[ F(x) = \frac{1}{1+(x/\alpha)^{-\beta}} \]

As an exercise, you can derive the density in the following form

\[ f(x) = \frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1+(x/\alpha)^\beta]^2}, \ \ \text{ for } \ \ x > 0. \]

After some algebra, we can find the \(k\)th moment

\[ \mu_k = E[X^k] = \alpha^k B\left(1+\frac{k}{\beta}, 1 - \frac{k}{\beta} \right). \]

This assignment will focus on finding MME of parameters \(\alpha\) and \(\beta\) based on a real-world application data set.

Question 1: Derive the log-logistic density function

Given the CDF of the two-parameter log-logistic distribution

\[ F(x) = \frac{1}{1+(x/\alpha)^{-\beta}}. \]

The distribution function of the log-logistic function is given above. Differentiation can be used to find the density function, f(x), of the log-logistic distribution.

f(x) = \(\frac{d}{dx}\) F(x) = \(\frac{d}{dx}\) \(\frac{1}{1+(x/\alpha)^{-\beta}}\)

This entire distribution function F(x) can be multiplied by \(\frac{(x/\alpha)^{\beta}}{(x/\alpha)^{\beta}}\) to remove the negative exponent in the denominator.

F(x) = \(\frac{1}{1+(x/\alpha)^{-\beta}}\) * \(\frac{(x/\alpha)^{\beta}}{(x/\alpha)^{\beta}}\) = \(\frac{(x/\alpha)^{\beta}}{(x/\alpha)^{\beta} + [ (x/\alpha)^{\beta}*(x/\alpha)^{-\beta}]}\) = \(\frac{(x/\alpha)^{\beta}}{(x/\alpha)^{\beta} + 1}\) = \(\frac{(x/\alpha)^{\beta}}{1 +(x/\alpha)^{\beta}}\)

So, F(x) = \(\frac{(x/\alpha)^{\beta}}{1 +(x/\alpha)^{\beta}}\) and f(x) = \(\frac{d}{dx}\) \(\frac{(x/\alpha)^{\beta}}{1 +(x/\alpha)^{\beta}}\)

In order to find the density function, differentiation with the quotient rule will be conducted. f(x) = \(\frac{u'(x)v(x) - u(x)v'(x)}{[v(x)]^2}\)

For this function, the parts of the quotient rule are as follows.

u(x) = \((x/\alpha)^{\beta}\) and u’(x) = \((\beta/\alpha) *(x/\alpha)^{\beta-1}\)

v(x) = \(1 +(x/\alpha)^{\beta}\) and v’(x) = \((\beta/\alpha) *(x/\alpha)^{\beta-1}\)

So, now the quotient rule formula will be used to find the density function f(x) for this log-logistic distribution.

f(x) = \(\frac{[(\beta/\alpha)(x/\alpha)^{\beta-1}(1 +(x/\alpha)^{\beta})] - [(x/\alpha)^{\beta}(\beta/\alpha) *(x/\alpha)^{\beta-1}]}{[1 +(x/\alpha)^{\beta}]^2}\) = \(\frac{[(\beta/\alpha)(x/\alpha)^{\beta-1} +(x/\alpha)^{\beta}(\beta/\alpha)(x/\alpha)^{\beta-1}] - [(x/\alpha)^{\beta}(\beta/\alpha) *(x/\alpha)^{\beta-1}]}{[1 +(x/\alpha)^{\beta}]^2}\)

= \(\frac{(\beta/\alpha)(x/\alpha)^{\beta-1} +(\beta/\alpha)(x/\alpha)^{2\beta-1} - (\beta/\alpha) *(x/\alpha)^{2\beta-1}}{[1 +(x/\alpha)^{\beta}]^2}\) = \(\frac{(\beta/\alpha)(x/\alpha)^{\beta-1} +(\beta/\alpha)(x/\alpha)^{2\beta-1} - (\beta/\alpha) *(x/\alpha)^{2\beta-1}}{[1 +(x/\alpha)^{\beta}]^2}\) = \(\frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1 +(x/\alpha)^{\beta}]^2}\), for \(x > 0\)

So, the density function, f(x), for the log-logistic distribution is given as:

f(x) = \(\frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1 +(x/\alpha)^{\beta}]^2}\) , for \(x > 0\)

Question 2: Distribution of Recovery Time from A Surgery

Time to recovery (in days) after a specific knee surgery procedure. This follows a typical log-logistic pattern in medical survival/recovery analysis:

8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56, 
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13, 
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33, 
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65

Based on the above data to perform the following analysis.

1. Using method of moment estimation to estimate \(\alpha\) and \(\beta\), denoted by \(\hat{\alpha}\) and \(\hat{\beta}\), respectively.

To begin, a data set of the above values will be created. These values represent the time to recovery, in days, after a knee surgery procedure. This data set will be called ‘recovery’.

recovery <- c(8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56, 27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13, 35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33, 44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65)

Next, the method of moments estimation will be used to estimate \(\alpha\) and \(\beta\) for the recovery data set. These estimated values will be denoted as \(\hat{\alpha}\) and \(\hat{\beta}\).

m1 <- mean(recovery)
m2 <- mean(recovery^2)

mom_equation <- function(beta) {
  term1 <- beta(1 + 2/beta, 1 - 2/beta)
  term2 <- beta(1 + 1/beta, 1 - 1/beta)^2
  return(m2/m1^2 - term1/term2)}

mom_beta <- uniroot(mom_equation, lower = 2, upper = 50)$root

mom_alpha <- m1 / beta(1 + 1/mom_beta, 1 - 1/mom_beta)

cat("MoM estimates:\n",
    "   alpha_hat =", round(mom_alpha, 4),
    ", beta_hat =", round(mom_beta, 4), "\n")

MoM estimates:
    alpha_hat = 32.6543 , beta_hat = 6.0062 

The estimates found by the method of moment estimation are as follows: \(\hat{\alpha}\) = 32.6543 and \(\hat{\beta}\) = 6.0062.

2. Since the moment estimates \(\hat{\alpha}\) and \(\hat{\beta}\) are random, construct bootstrap sampling distributions for each. To visualize these distributions, plot separate bootstrap histograms for \(\hat{\alpha}\) and \(\hat{\beta}\). Then, overlay a smooth density curve on each histogram using Gaussian kernel density estimation. Finally, describe the patterns of these density curves.

I will begin by generating the bootstrap estimates to construct the bootstrap sampling distribution.

set.seed(123)
B <- 5000

alpha_boot <- numeric(B)
beta_boot <- numeric(B)
n <- length(recovery)

for (i in 1:B) {
  boot_sample <- sample(recovery, n, replace = TRUE)
  
  # The first and second moments
  m1_boot <- mean(boot_sample)
  m2_boot <- mean(boot_sample^2)
  
  # The method of moments equation
  mom_equation <- function(beta) {
    term1 <- beta(1 + 2/beta, 1 - 2/beta)
    term2 <- beta(1 + 1/beta, 1 - 1/beta)^2
    return(m2_boot/m1_boot^2 - term1/term2)
  }
  
  # Solving for beta hat
  beta_hat <- tryCatch(
    uniroot(mom_equation, lower = 2, upper = 50)$root,
    error = function(e) NA )
  
  # Solving for alpha hat
    alpha_hat <- m1_boot / beta(1 + 1/beta_hat, 1 - 1/beta_hat)
    alpha_boot[i] <- alpha_hat
    beta_boot[i] <- beta_hat
}

# I will remove missing values
alpha_boot <- alpha_boot[!is.na(alpha_boot)]
beta_boot <- beta_boot[!is.na(beta_boot)]

Now, I will plot the bootstrap distributions of \(\hat{\alpha}\) and \(\hat{\beta}\). For each plot, I will overlay a smooth density curve on each histogram using Gaussian kernel density estimation.

par(mfrow = c(1,2))

# Alpha Hat
hist(alpha_boot, probability = TRUE, breaks = 25,
     main = "Bootstrap for Alpha Hat",
     xlab = "Alpha Hat")
lines(density(alpha_boot), col = "purple", lwd = 2)

# Beta Hat
hist(beta_boot, probability = TRUE, breaks = 25,
     main = "Bootstrap for beta Hat",
     xlab = "Beta Hat")
lines(density(beta_boot), col = "purple", lwd = 2)

The bootstrap distribution of \(\hat{\alpha}\) appears to be symmetric and unimodal, with an approximately Normal distribution. The mean of \(\hat{\alpha}\) is roughly between 32 and 33, which matches up perfectly with the estimate of \(\hat{\alpha}\) from the previous part which was found to be 32.6543. The bootstrap distribution of \(\hat{\beta}\) is unimodal and appears to be ever so slightly skewed to the right. The mean of the \(\hat{\beta}\) appears to be around 6, which matches up perfectly with the estimate of \(\hat{\beta}\) = 6.0062 from the previous part, which was found to be 6.0062.