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

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

Given the CDF of log-logistic distribution: \[ F(x)=\frac{1}{1+(x/\alpha)^{-\beta}},\quad x>0. \]

Rewrite it as: \[ F(x)=\frac{(x/\alpha)^\beta}{1+(x/\alpha)^\beta}. \]

Let \(t=(x/\alpha)^\beta\). Then \(F(x)=t/(1+t)\), so \[ f(x)=F'(x)=\frac{t'}{(1+t)^2}. \]

Since \(t'=\frac{\beta}{\alpha}(x/\alpha)^{\beta-1}\), we obtain

\[ f(x)=\frac{\beta}{\alpha}\frac{(x/\alpha)^{\beta-1}}{\left[1+(x/\alpha)^\beta\right]^2},\quad 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.

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)

# sample moments
m1 <- mean(x)
m2 <- mean(x^2)

# target ratio m2/m1^2
target <- m2 / (m1^2)

# MoM equation after eliminating alpha
g_beta <- function(beta){
  A1 <- beta(1 + 1/beta, 1 - 1/beta)
  A2 <- beta(1 + 2/beta, 1 - 2/beta)
  (A2 / (A1^2)) - target
}

# beta must be > 2 (so second moment exists)
beta_hat <- uniroot(g_beta, interval = c(2.01, 50))$root

# plug-in alpha_hat
A1_hat <- beta(1 + 1/beta_hat, 1 - 1/beta_hat)
alpha_hat <- m1 / A1_hat

# # one lean output box: MoM point estimates
data.frame(alpha_hat = alpha_hat, beta_hat = beta_hat)

  alpha_hat beta_hat
1   32.6543 6.006232

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.

set.seed(123)

B <- 1000

alpha_star <- NULL
beta_star  <- NULL

for(i in 1:B){
  
  xb <- sample(x, size = n, replace = TRUE)
  
  m1b <- mean(xb)
  m2b <- mean(xb^2)
  target_b <- m2b / (m1b^2)
  
  # MoM equation for this bootstrap sample
  g_beta_b <- function(beta){
    A1 <- beta(1 + 1/beta, 1 - 1/beta)
    A2 <- beta(1 + 2/beta, 1 - 2/beta)
    (A2 / (A1^2)) - target_b
  }
  
  # Solve for beta_hat* 
  out <- try(uniroot(g_beta_b, interval = c(2.01, 50))$root, silent = TRUE)
  
  if(inherits(out, "try-error")){
    beta_star[i]  <- NA
    alpha_star[i] <- NA
  } else {
    beta_star[i] <- out
    A1b <- beta(1 + 1/beta_star[i], 1 - 1/beta_star[i])
    alpha_star[i] <- m1b / A1b
  }
}

# remove NA's 
alpha_star <- alpha_star[!is.na(alpha_star)]
beta_star  <- beta_star[!is.na(beta_star)]

# one lean output box: bootstrap summary
data.frame(
  B_used = length(alpha_star),
  alpha_boot_mean = mean(alpha_star),
  alpha_boot_sd = sd(alpha_star),
  beta_boot_mean = mean(beta_star),
  beta_boot_sd = sd(beta_star)
)

  B_used alpha_boot_mean alpha_boot_sd beta_boot_mean beta_boot_sd
1   1000        32.65709      1.644027       6.134973    0.5973591

par(mfrow = c(1,2), cex.main = 0.85)

hist(alpha_star, probability = TRUE, breaks = 20,
     main = "Bootstrap of alpha_hat",
     xlab = expression(hat(alpha)))
lines(density(alpha_star), col = "blue", lwd = 2)

hist(beta_star, probability = TRUE, breaks = 20,
     main = "Bootstrap of beta_hat",
     xlab = expression(hat(beta)))
lines(density(beta_star), col = "blue", lwd = 2)

par(mfrow = c(1,1))

Answer: The bootstrap distribution of \(\hat{\alpha}\) is unimodal and roughly symmetric, with most values concentrated around 32–33, suggesting that the MoM estimate of \(\alpha\) is fairly stable. In contrast, the bootstrap distribution of \(\hat{\beta}\) is unimodal but clearly right-skewed: it peaks near 6 and has a longer right tail extending toward about 8–9. This indicates larger sampling variability for \(\hat{\beta}\), which is reasonable because \(\hat{\beta}\) is obtained by numerically solving a nonlinear moment equation.