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}}. \]

# Install packages

library(Deriv)

# Define the CDF
F_fun <- function(x, alpha, beta){
  1/(1 + (x/alpha)^(-beta))
}

# Differentiate to x
f_fun <- Deriv(F_fun, "x")

f_fun

function (x, alpha, beta) 
{
    .e1 <- x/alpha
    beta/(alpha * (1 + 1/.e1^beta)^2 * .e1^(1 + beta))
}

# Define the log-logistic PDF
f_loglogistic <- function(x, alpha, beta){
  
  # Implements simplified formula for the density
  (beta/alpha) * (x/alpha)^(beta - 1) /
    (1 + (x/alpha)^beta)^2
}

# Parameter choices for illustration
alpha <- 30
beta  <- 3

# Create a sequence of x values
x <- seq(0.1, 80, length = 400)

# Plot the density curve
plot(x, f_loglogistic(x, alpha, beta),
     type = "l",
     lwd = 2,
     main = "Log-Logistic Density Function Visualization",
     xlab = "x",
     ylab = "f(x)")

Starting from: \[ F(x) = \frac{1}{1+(x/\alpha)^{-\beta}}. \] Differentiating with respect to x gives: \[ f(x) = \frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1+(x/\alpha)^\beta]^2} \]

which is the log-logistic probability density function.

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.

# Data
x <- 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
)

n <- length(x)

# Compute sample moments
m <- mean(x)
s2 <- var(x)

m

[1] 34.1922

s2

[1] 122.1819

# Beta
moment_eq <- function(beta){
  
  A <- (pi/beta)/sin(pi/beta)
  B <- (2*pi/beta)/sin(2*pi/beta)
  
  alpha <- m / A
  
  theoretical_var <- alpha^2 * (B - A^2)
  
  theoretical_var - s2
}

beta_hat <- uniroot(moment_eq, interval = c(2.01, 20))$root

beta_hat

[1] 5.95252

# Alpha 
alpha_hat <- m * sin(pi/beta_hat)/(pi/beta_hat)

alpha_hat

[1] 32.62681

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.

# Bootstrap Sampling Distributions
set.seed(1)

B <- 2000

alpha_boot <- numeric(B)
beta_boot  <- numeric(B)

for(i in 1:B){
  
  xb <- sample(x, n, replace = TRUE)
  
  m_b <- mean(xb)
  s2_b <- var(xb)
  
  moment_eq_boot <- function(beta){
    
    A <- (pi/beta)/sin(pi/beta)
    Bterm <- (2*pi/beta)/sin(2*pi/beta)
    
    alpha <- m_b / A
    
    theoretical_var <- alpha^2 * (Bterm - A^2)
    
    theoretical_var - s2_b
  }
  
  beta_hat_b <- tryCatch(
    uniroot(moment_eq_boot, interval = c(2.01, 20))$root,
    error = function(e) NA
  )
  
  if(!is.na(beta_hat_b)){
    alpha_hat_b <- m_b * sin(pi/beta_hat_b)/(pi/beta_hat_b)
  } else {
    alpha_hat_b <- NA
  }
  
  alpha_boot[i] <- alpha_hat_b
  beta_boot[i]  <- beta_hat_b
}

alpha_boot <- na.omit(alpha_boot)
beta_boot  <- na.omit(beta_boot)

# Bootstrap Histogram for Alpha
hist(alpha_boot,
     probability = TRUE,
     main = expression("Bootstrap Distribution of " * hat(alpha)),
     xlab = expression(hat(alpha)))

lines(density(alpha_boot), lwd = 2)

# Bootstrap Histogram for Beta
hist(beta_boot,
     probability = TRUE,
     main = expression("Bootstrap Distribution of " * hat(beta)),
     xlab = expression(hat(beta)))

lines(density(beta_boot), lwd = 2)

The bootstrap distribution of the scale estimator of \(\alpha\) and \(\beta\) are both symmetric and unimodal and shows that bootstrap resampling provides a accurate approximation of the estimators’ sampling variability.