Assignment 9: Likelihood-based Chi-Square Tests

Assignment Objectives

• Enhance understanding the procedure of likelihood-based chi-square hypothesis testing .

• Implement the procedures for detecting overfitting/underfitting issues in practical applications.

Policies of Using AI Tools

Policy on AI Tool Use: Please 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.

Testing Overfitting/Underfitting

In Machine Learning and Statistics, overfitting occurs when a model is too complex and learns noise, leading to poor performance on new data, while underfitting happens when a model is too simple to capture important patterns, resulting in high errors overall; both issues are explained by the Bias–Variance Tradeoff and can cause unreliable predictions in real-world applications.

The probability density function (PDF) of the Weibull distribution is:

\[ f(t; \lambda, \beta) = \frac{\beta}{\lambda} \left( \frac{t}{\lambda} \right)^{\beta-1} \exp\left[ -\left( \frac{t}{\lambda} \right)^\beta \right], \quad t \ge 0 \] where \(\lambda > 0\) is the scale parameter (characteristic life) and \(\beta > 0\) is the shape parameter.

When \(\beta = 1\), the Weibull PDF simplifies to the exponential PDF:

\[ f(t; \lambda) = \frac{1}{\lambda} \exp\left( -\frac{t}{\lambda} \right) \] with constant hazard rate \(h(t) = 1/\lambda\).

This assignment focuses on performing a hypothesis test for the shape parameter (\(\beta\)) of the Weibull distribution within a reliability mode

\[\begin{align} H_0&: \beta = 1 \quad \text{(Exponential model, simpler)} \\ H_1&: \beta \neq 1 \quad \text{(Weibull model, more complex)} \end{align}\]

Question: Reliability Application

A mid-sized manufacturing company producing industrial conveyor systems began experiencing unexpected downtime in one of its distribution facilities, prompting concern about the reliability of a newly sourced batch of ball bearings used in the motor assemblies. These bearings, supplied by a vendor adopting cost-saving production methods, were installed across multiple units operating under continuous load conditions. After several months, maintenance logs revealed a pattern of increasing failures, with components lasting anywhere from a few dozen to over 150 hours before breakdown. To investigate, the engineering team collected time-to-failure data from 50 identical bearings and conducted a Weibull analysis within the framework of Reliability Engineering. The 50 time-to-failure (survival time) are:

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

# Store the 50 failure times as a numeric vector for analysis
time <- 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
)

This assignment focuses on hypothesis \(H_0: \beta = 1\) (exponential) against \(H_1: \beta \neq 1\) (Weibull). This framework detects overfitting (fitting a Weibull when exponential is true) and underfitting (fitting exponential when Weibull with \(\beta \neq 1\) is true).

a). Find the MLE of the Weibull parameters \(\lambda\) (scale) and \(\beta\) (shape), denoted by \(\hat{\lambda}\) and \(\hat{\beta}\), respectively, using the optim() procedure. [Hint: You should provide explicit expressions for the log-likelihood and gradient functions of the Weibull distribution parameters.]

# Define the Weibull log-likelihood for (lambda, beta)
weibull_loglik <- function(par, x) {
  lambda <- par[1]
  beta <- par[2]
  
  # Restrict parameters to positive values
  if (lambda <= 0 || beta <= 0) return(-Inf)
  
  n <- length(x)
  
  n * log(beta) -
    n * beta * log(lambda) +
    (beta - 1) * sum(log(x)) -
    sum((x / lambda)^beta)
}

# Define the gradient for faster optimization
weibull_grad <- function(par, x) {
  lambda <- par[1]
  beta <- par[2]
  n <- length(x)
  
  # Repeated term in derivatives
  s <- (x / lambda)^beta
  
  d_lambda <- (beta / lambda) * (sum(s) - n)
  d_beta <- n / beta - n * log(lambda) + sum(log(x)) -
    sum(s * log(x / lambda))
  
  c(d_lambda, d_beta)
}

# Use optim() to maximize the log-likelihood numerically
fit_weibull <- optim(
  par = c(mean(time), 1.5),   # starting values
  fn = function(par, x) -weibull_loglik(par, x),  # minimize negative log-likelihood
  gr = function(par, x) -weibull_grad(par, x),    # negative gradient
  x = time,
  method = "BFGS",
  hessian = TRUE
)

# Extract MLE estimates
fit_weibull$par

[1] 98.992745  2.205676

Using the optim() procedure to maximize the Weibull log-likelihood, the maximum likelihood estimates are \(\hat{\lambda} = 98.992745\) and \(\hat{\beta} = 2.205676\).

b). Find the MLE of the exponential parameter \(\lambda\) (scale), denoted by \(\hat{\lambda}\), using any procedure. [Hint: You should provide explicit expressions for the log-likelihood and gradient functions of the exponential distribution parameters.]

For the exponential model, the log-likelihood is \(\ell(\lambda) = -n \log \lambda - \frac{1}{\lambda} \sum_{i=1}^n t_i\).

The gradient is \(\frac{d\ell}{d\lambda} = -\frac{n}{\lambda} + \frac{1}{\lambda^2} \sum_{i=1}^n t_i\).

Setting the derivative equal to zero gives \(\hat{\lambda} = \bar{t}\).

# Define the exponential log-likelihood for lambda
exp_loglik <- function(lambda, x) {
  n <- length(x)
  -n * log(lambda) - sum(x) / lambda
}

# Define the gradient of the exponential log-likelihood
exp_grad <- function(lambda, x) {
  n <- length(x)
  -n / lambda + sum(x) / lambda^2
}

# Compute the exponential MLE using the sample mean
lambda_hat_exp <- mean(time)

lambda_hat_exp

[1] 87.61

Using the sample mean, the maximum likelihood estimate is \(\hat{\lambda} = 87.61\).

c). Perform the likelihood ratio \(\chi^2\) test on \(\beta = 1\). What is the p-value? [Hint: Use the results in a) and b).]

The likelihood ratio statistic is \(LR = 2\left[\ell(\hat{\lambda}, \hat{\beta}) - \ell(\hat{\lambda}_0, 1)\right]\).

# Compute log-likelihood for Weibull (full model)
loglik_weibull <- weibull_loglik(fit_weibull$par, time)

# Compute log-likelihood for exponential (restricted model, beta = 1)
lambda_hat_exp <- mean(time)
n <- length(time)

loglik_exp <- -n * log(lambda_hat_exp) - sum(time) / lambda_hat_exp

# Compute likelihood ratio statistic
LR <- 2 * (loglik_weibull - loglik_exp)

# Compute p-value from chi-square distribution with 1 df
p_value_LR <- pchisq(LR, df = 1, lower.tail = FALSE)

LR

[1] 34.44994

p_value_LR

[1] 4.373569e-09

The likelihood ratio statistic is \(LR = 34.44994\). The corresponding p-value is approximately \(4.373569 \times 10^{-9}\).

d). Perform the Wald \(\chi^2\) on the same hypothesis \(\beta = 1\). What is the p-value? [Hint: You need to find the observed Fisher information matrix (i.e., the negative Hessian matrix from optim()) based on the Weibull distribution. The inverse of the negative observed Hessian matrix is the covariance matrix.]

The Wald statistic for testing \(H_0:\beta = 1\) is \(W = \dfrac{(\hat{\beta} - 1)^2}{\widehat{\mathrm{Var}}(\hat{\beta})}\).

# Extract the Hessian from the Weibull fit
info_weibull <- fit_weibull$hessian

# Invert the Hessian to get the covariance matrix
cov_weibull <- solve(info_weibull)

# Extract the Weibull estimate and its variance
beta_hat <- fit_weibull$par[2]
var_beta_hat <- cov_weibull[2, 2]

# Compute the Wald chi-square statistic
W <- (beta_hat - 1)^2 / var_beta_hat

# Compute the p-value using a chi-square distribution with 1 df
p_value_W <- pchisq(W, df = 1, lower.tail = FALSE)

W

[1] 23.10267

p_value_W

[1] 1.535777e-06

The Wald statistic is \(W = 23.10267\). The corresponding p-value is \(1.535777 \times 10^{-6}\).

e). Write a summary of the above analyses to address the following:

• Whether the two tests generated the same results.

• Which model is recommended for the data.

• Draw the density curve based on the MLE(s) of the parameter(s) and describe the distribution of the time-to-failure.

Both the likelihood ratio test and the Wald test generate the same result. The likelihood ratio test gave a p-value of \(4.373569 \times 10^{-9}\), and the Wald test gave a p-value of \(1.535777 \times 10^{-6}\). Both p-values are extremely small and are statistically significant. We have enough evidence to reject the null hypothesis \(H_0: \beta = 1\).

The null hypothesis represents the exponential model, so rejecting \(H_0\) means that the exponential model does not fit this data well. Because of this, the Weibull model is a better choice for modeling the time-to-failure data.

# Draw the fitted Weibull density using the MLEs from part (a)
curve(
  dweibull(x, shape = fit_weibull$par[2], scale = fit_weibull$par[1]),
  from = 0, to = 200,
  lwd = 2,
  xlab = "Time to failure",
  ylab = "Density",
  main = "Fitted Weibull Density"
)

Using the MLEs from part (a), the fitted Weibull density is unimodal and right-skewed. The estimated Weibull shape parameter is \(\hat{\beta} = 2.205676\), which is greater than 1. This suggests that the failure rate is not constant, but instead increases as time goes on. The fitted Weibull density shows that larger time values are associated with a higher concentration of failures. Based on our observation, this is consistent with wear-out failure in reliability analysis. Components tend to weaken over time, so the longer they are in use, the more likely they are to experience failure.