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

Answer to Question 1

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

Answer to Question 2 Part A:

Formula for k-th moment: \[ \mu_k = E[X^k] = \alpha^k B\left(1+\frac{k}{\beta}, 1 - \frac{k}{\beta} \right)=\alpha^k\cdot\frac{\Gamma(1+k/\beta)\cdot\Gamma(1-k/\beta)}{\Gamma(1+k/\beta+1-k/\beta)}=\alpha^k\cdot\Gamma(1+k/\beta)\cdot\Gamma(1-k/\beta) \]

First and second population moments:

\[ \mu_1=E[X]=\alpha\cdot\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta)\\\\ \mu_2=E[X^2]=\alpha^2\cdot\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta) \]

First and second sample moments:

\[ m_1=\frac{1}{n}\sum_{i=1}^{n} X_i\\\\ m_2=\frac{1}{n}\sum_{i=1}^{n} X_i^2 \]

Method of Moments Equations:

\[ m_1=\alpha\cdot\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta)\\\\ m_2=\alpha^2\cdot\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta) \]

Finding equations to solve numerically for \(\alpha\) and \(\beta\):

\[ m_1^2=\alpha^2\cdot(\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta))^2\\\\ \frac{m_2}{m_1^2}=\frac{\alpha^2\cdot\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta)}{\alpha^2\cdot(\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta))^2}=\frac{\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta)}{(\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta))^2} \]

Let \(\hat{\beta}\) be the solution to:

\[ \frac{m_2}{m_1^2}=\frac{\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta)}{(\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta))^2} \]

Let \(\hat{\alpha}\) be obtained from:

\[ \hat{\alpha}=\frac{m_1}{\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta)} \]

The parameter \(\beta\) can be solved numerically from the following equation:

\[ g(\beta)=\frac{\Gamma(1+2/\beta)\cdot\Gamma(1-2/\beta)}{(\Gamma(1+1/\beta)\cdot\Gamma(1-1/\beta))^2}-\frac{m_2}{m_1^2}=0 \] Plotting \(g(\beta)\) to identify the interval where the root is located:

#Setting the sample of the recovery time in days (data is above) equal to x.
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)

#Finding 1st and 2nd moments of sample.
m1 <- mean(x)
m2 <- mean(x^2)

#Allows us to use g(beta) function from above in R (here using b for simplicity)
gb <- function(b){((gamma(1+2/b))*(gamma(1-2/b))/(((gamma(1+1/b))*(gamma(1-1/b)))^2)) - m2/(m1^2)}

#Creating data frame with x values between 4 and 10 and y values calculated from gb based on those x values.
x_vals <- seq(4,10,length.out=50)
dframe <- data.frame(x=x_vals, y=gb(x_vals))

#Creates plot of gb values from dframe with plotted x and y axis.
gb.plot <- ggplot(dframe, aes(x = x, y = y)) +
  geom_line(color="steelblue", size= 1) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  labs(title = "The curve of function g(beta)",
       x = "beta",
       y = "g(beta)") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.margin = margin(t = 35, r = 20, l = 30, unit = "pt"))
ggplotly(gb.plot)

The curve above indicates that the root lies in (5,7.5). We can use this to find beta and alpha:

library(stats)
#Using uniroot() to find beta
beta <- uniroot(gb, interval = c(5, 7.5))$root
#Using beta to find alpha from equation determing above.
alpha <- m1/((gamma(1+1/beta))*(gamma(1-1/beta)))
pander(cbind(alpha = alpha, beta = beta))

alpha	beta
32.65	6.006

Estimating \(\alpha\) and \(\beta\), denoted by \(\hat{\alpha}\) and \(\hat{\beta}\), respectively, we get the result \(\hat{\alpha}=32.65\) and \(\hat{\beta}=6.006\)

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}\). hen, overlay a smooth density curve on each histogram using Gaussian kernel density estimation. Finally, describe the patterns of these density curves.

Answer to Question 2 Part B:

set.seed(123) ##Used to set seed so that result can be replicated

optim_mom <- function(data){ #Creates function that can be used later
  m1 <- mean(data)
  m2 <- mean(data^2)
  target <- m2/(m1^2) #Used to get part of objective function
  
  
  objective <- function(k){
    part1 <- gamma(1+2/k)*gamma(1-2/k)
    part2 <- gamma(1+1/k)*gamma(1-1/k)
    (part1/(part2^2)-target)^2    #Objective function to minimize                   
  }
  #Optimizes for k using Brent's Method to get minimized value
  opt <- optimize(objective, interval=c(0.1,10))
  beta_hat <- opt$minimum
  alpha_hat <- m1/((gamma(1+1/beta_hat))*(gamma(1-1/beta_hat)))
  
  return(list(alpha = alpha_hat, beta = beta_hat))
}


n <- length(x) #Sets sample size for bootstrap equal to length of original sample
B <- 1000 #How many times we should bootstrap sample

bootstrap.beta <- numeric(B) #Resets 'bootstrap.beta' 
bootstrap.alpha <- numeric(B) #Resets 'bootstrap.alpha'



for (i in 1:B){
  boot.sample <- sample(x, size=n, replace=TRUE)
  result <- optim_mom(boot.sample) #Uses function created above and stores results in 'result'
  bootstrap.beta[i] <- result$beta
  bootstrap.alpha[i] <- result$alpha
  #Used to randomly sample (with replacement) from 'x' 'B' times to get sample beta-hats and alpha-hats
}

hist(bootstrap.alpha, prob=TRUE, #Creates histogram of bootstrap alpha-hats
     main = "Bootstrap Distribution of Alpha-Hat",
     xlab = "Alpha-Hats")
lines(density(bootstrap.alpha, kernel = "gaussian", bw=1), col="steelblue", lwd=2) #Creates density curve usin Guassian kernel

hist(bootstrap.beta, prob=TRUE, #Creates histogram of bootstrap beta-hats
     main = "Bootstrap Distribution of Beta-Hats",
     xlab = "Beta-Hats")
lines(density(bootstrap.beta, kernel = "gaussian", bw=0.2), col="steelblue", lwd=2) #Creates density curve usin Guassian kernel

The histogram for the bootstrap \(\hat{\alpha}\)’s is centered somewhere between 32 and 34 (which makes sense given that the \(\hat{\alpha}\) for the sample was 32.65). The density curve for this histogram appears to follow a symmetric distribution. Alternatively, the histogram for the bootstrap \(\hat{\beta}\) is centered around 6 (which again makes sense given that the \(\hat{\beta}\) for the sample was 6.006). The density curve for this histogram appears to be skewed slightly to the right.