STA 506 HOMEWORK 10

1 Question A MLE of Weibull Parameters

The Weibull probability density function is:

\[ f(t;\lambda,\beta) = \frac{\beta}{\lambda} \left(\frac{t}{\lambda}\right)^{\beta-1} \exp\left[-\left(\frac{t}{\lambda}\right)^\beta\right] \]

For an independent sample \(t_1, \dots, t_n\), the log-likelihood function is:

\[ \ell(\lambda,\beta) = n \ln \beta - n \beta \ln \lambda + (\beta - 1)\sum_{i=1}^n \ln t_i - \sum_{i=1}^n \left(\frac{t_i}{\lambda}\right)^\beta \]

The gradient components are:

\[ \frac{\partial \ell}{\partial \lambda} = -\frac{n\beta}{\lambda} + \frac{\beta}{\lambda} \sum_{i=1}^n \left(\frac{t_i}{\lambda}\right)^\beta \]

\[ \frac{\partial \ell}{\partial \beta} = \frac{n}{\beta} - n \ln \lambda + \sum_{i=1}^n \ln t_i - \sum_{i=1}^n \left(\frac{t_i}{\lambda}\right)^\beta \ln\left(\frac{t_i}{\lambda}\right) \]

Since no closed-form solution exists, numerical optimization is used.

t <- c(12.4, 18.7, 25.3, 30.1, 33.5, 35.2, 38.9, 40.3, 42.7, 45.1,
       47.6, 49.8, 52.4, 55.0, 57.3, 60.2, 62.8, 65.1, 67.9, 70.5,
       72.3, 75.6, 78.2, 80.9, 83.4, 85.7, 88.1, 90.6, 93.2, 95.8,
       98.4, 101.0, 104.5, 107.3, 110.6, 113.2, 116.8, 120.1, 123.7,
       127.4, 130.9, 134.5, 138.2, 142.0, 146.3, 150.7, 155.2, 160.8,
       168.4, 175.9)

n <- length(t)

# Log-likelihood function (negative for minimization)
loglik_weibull <- function(par) {
  lambda <- par[1]
  beta <- par[2]
  
  if(lambda <= 0 || beta <= 0) return(Inf)
  
  val <- n*log(beta) - n*beta*log(lambda) +
         (beta - 1)*sum(log(t)) -
         sum((t/lambda)^beta)
  
  return(-val)  # negative for minimization
}

# Gradient
grad_weibull <- function(par) {
  lambda <- par[1]
  beta <- par[2]
  
  d_lambda <- (-n*beta/lambda) +
              (beta/lambda)*sum((t/lambda)^beta)
  
  d_beta <- n/beta - n*log(lambda) +
            sum(log(t)) -
            sum((t/lambda)^beta * log(t/lambda))
  
  return(-c(d_lambda, d_beta))  # negative for minimization
}

# Initial guesses
init <- c(mean(t), 1)

# Optimization
fit <- optim(init, loglik_weibull, grad_weibull, method = "BFGS", hessian = TRUE)

lambda_hat <- fit$par[1]
beta_hat <- fit$par[2]

lambda_hat

[1] 98.95336

beta_hat

[1] 2.205325

2 Question B MLE of Exponential Parameter \(\lambda\)

When \(\beta = 1\), the Weibull distribution reduces to the exponential distribution:

\[ f(t;\lambda) = \frac{1}{\lambda} \exp\left(-\frac{t}{\lambda}\right) \]

For an independent sample \(t_1, \dots, t_n\), the log-likelihood function is:

\[ \ell(\lambda) = -n \ln \lambda - \frac{1}{\lambda}\sum_{i=1}^n t_i \]

Taking the derivative with respect to \(\lambda\) and setting it equal to zero:

\[ \frac{d\ell}{d\lambda} = -\frac{n}{\lambda} + \frac{1}{\lambda^2}\sum_{i=1}^n t_i = 0 \]

Solving for \(\lambda\), we obtain:

\[ \hat{\lambda} = \bar{t} \]

# Data (reuse)
t <- c(12.4, 18.7, 25.3, 30.1, 33.5, 35.2, 38.9, 40.3, 42.7, 45.1,
       47.6, 49.8, 52.4, 55.0, 57.3, 60.2, 62.8, 65.1, 67.9, 70.5,
       72.3, 75.6, 78.2, 80.9, 83.4, 85.7, 88.1, 90.6, 93.2, 95.8,
       98.4, 101.0, 104.5, 107.3, 110.6, 113.2, 116.8, 120.1, 123.7,
       127.4, 130.9, 134.5, 138.2, 142.0, 146.3, 150.7, 155.2, 160.8,
       168.4, 175.9)

# MLE of lambda for exponential distribution
lambda_hat_exp <- mean(t)

lambda_hat_exp

[1] 87.61

3 Question C Likelihood Ratio \(\chi^2\) Test for \(H_0: \beta = 1\)

We test:

\[ H_0: \beta = 1 \quad \text{vs} \quad H_1: \beta \ne 1 \]

The likelihood ratio statistic is defined as:

\[ LR = 2\big[\ell(\hat{\lambda}, \hat{\beta}) - \ell(\hat{\lambda}_{H_0}, 1)\big] \]

where: - \(\ell(\hat{\lambda}, \hat{\beta})\) is the log-likelihood under the full Weibull model - \(\ell(\hat{\lambda}_{H_0}, 1)\) is the log-likelihood under the exponential model

By Wilks’ Theorem:

\[ LR \xrightarrow{d} \chi^2_1 \]

The p-value is computed as:

\[ p\text{-value} = P(\chi^2_1 \ge LR) \]

# Reuse results from (a) and (b)
# theta_hat values already computed:
lambda_hat_weibull <- lambda_hat
beta_hat <- beta_hat

lambda_hat_exp <- mean(t)

n <- length(t)

# Log-likelihood for Weibull (full model)
loglik_weibull_full <- function(lambda, beta) {
  n*log(beta) - n*beta*log(lambda) +
  (beta - 1)*sum(log(t)) -
  sum((t/lambda)^beta)
}

# Log-likelihood for exponential (beta = 1)
loglik_exp <- function(lambda) {
  -n*log(lambda) - sum(t)/lambda
}

# Compute log-likelihoods
l_full <- loglik_weibull_full(lambda_hat_weibull, beta_hat)
l_restricted <- loglik_exp(lambda_hat_exp)

# Likelihood ratio statistic
LR <- 2 * (l_full - l_restricted)

# p-value
p_value <- 1 - pchisq(LR, df = 1)

p_value

[1] 4.373757e-09

4 Question D Wald \(\chi^2\) Test for \(H_0: \beta = 1\)

We test:

\[ H_0: \beta = 1 \quad \text{vs} \quad H_1: \beta \ne 1 \]

The Wald test statistic is given by:

\[ W = \frac{(\hat{\beta} - 1)^2}{\text{Var}(\hat{\beta})} \]

Under \(H_0\), the statistic follows:

\[ W \xrightarrow{d} \chi^2_1 \]

The variance of \(\hat{\beta}\) is obtained from the covariance matrix:

\[ \text{Cov}(\hat{\lambda}, \hat{\beta}) = H^{-1} \]

where \(H\) is the Hessian matrix obtained from the optimization procedure.

# Extract Hessian from optimization
H <- fit$hessian

# Covariance matrix (inverse of Hessian)
cov_matrix <- solve(H)

# Variance of beta_hat (second diagonal element)
var_beta <- cov_matrix[2,2]

# Wald statistic
W <- (beta_hat - 1)^2 / var_beta

# p-value
p_value_wald <- 1 - pchisq(W, df = 1)

p_value_wald

[1] 1.54359e-06

5 Question E Summary, Model Recommendation, and Density Visualization

5.1 Comparison of Tests

Two hypothesis tests were conducted to evaluate:

\[ H_0: \beta = 1 \quad \text{vs} \quad H_1: \beta \ne 1 \]

The results are:

• Likelihood Ratio Test p-value: \(4.37 \times 10^{-9}\)
• Wald Test p-value: \(1.54 \times 10^{-6}\)

Both p-values are significantly smaller than 0.05, leading to rejection of \(H_0\) in both cases.

Although both tests lead to the same decision, the likelihood ratio test produces a smaller p-value, indicating stronger evidence against the null hypothesis. This is consistent with theoretical expectations, as the likelihood ratio test is based directly on the likelihood function, whereas the Wald test relies on an asymptotic normal approximation.

---

5.2 Model Recommendation

The null hypothesis \(H_0: \beta = 1\) corresponds to the exponential model, which assumes a constant hazard rate. Rejecting this hypothesis indicates that this assumption is not supported by the data.

Since the exponential distribution is a special case of the Weibull distribution obtained by fixing \(\beta = 1\), rejecting \(H_0\) implies that this restriction is invalid. Therefore, the more general Weibull model, which allows \(\beta \ne 1\), provides a better fit to the data.

Based on the estimated value of \(\hat{\beta}\), the failure rate is not constant over time, further supporting the use of the Weibull model.

---

5.3 Density Plot

To visualize the fitted distribution, we plot the Weibull density using the maximum likelihood estimates \(\hat{\lambda}\) and \(\hat{\beta}\).

# Create sequence
t_seq <- seq(min(t), max(t), length.out = 200)

# Weibull density
weibull_density <- (beta_hat / lambda_hat_weibull) *
                   (t_seq / lambda_hat_weibull)^(beta_hat - 1) *
                   exp(-(t_seq / lambda_hat_weibull)^beta_hat)

# Set proper y-limit
y_max <- max(weibull_density)

# Plot histogram with adjusted ylim
hist(t, probability = TRUE, breaks = 10,
     ylim = c(0, y_max * 1.1),   # add buffer so curve fits
     main = "Weibull Fit to Time-to-Failure Data",
     xlab = "Time")

# Add curve
lines(t_seq, weibull_density, lwd = 2)

5.4 Interpretation of the Distribution

The fitted Weibull density curve aligns with the observed data and deviates from the shape expected under an exponential distribution. Specifically, the distribution does not exhibit the constant hazard rate implied by \(\beta = 1\).

Instead, the estimated shape parameter \(\hat{\beta}\) indicates that the failure rate changes over time. This behavior is consistent with the rejection of the exponential model and supports the conclusion that the Weibull distribution provides a more appropriate model for the time-to-failure data.