Assignment 12: Bootstrap Likelihood Ratio Test (BLRT)

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

Steps of the BLRT

• Fit models under \(H_0\) and \(H_1\)} to the original data, compute \(\Lambda_{\text{obs}}\).

• Generate bootstrap samples under \(H_0\)}:

  • Estimate parameters under \(H_0\) from the original data.
  • Generate \(B\) datasets by sampling from the model under \(H_0\) (parametric bootstrap) or by resampling residuals/cases (nonparametric bootstrap; parametric is common for BLRT).

• For each bootstrap sample \(b = 1,\dots,B\):

  • Fit \(H_0\) and \(H_1\) models.
  • Compute \(\Lambda_b = -2[\ell_{0,b} - \ell_{1,b}]\).

• Approximate p-value:

\[ p = \frac{1}{B} \sum_{b=1}^B I(\Lambda_b \ge \Lambda_{\text{obs}}) \] (Often a small adjustment is made for stability: \((1 + \#\{\Lambda_b \ge \Lambda_{\text{obs}}\})/(B+1)\)).

Question: Reliability Application

A wind energy company monitors the reliability of gearboxes in 75 identical wind turbines located in a coastal wind farm. The gearbox is a critical component; its failure often leads to costly downtime and repairs. Previous studies suggest that the hazard rate (failure risk) may increase over time due to mechanical wear (fatigue, pitting, bearing degradation). Engineers want to test whether the failure time distribution follows an exponential model (constant hazard, random failures) or a Weibull model with shape parameter \(k>1\) (increasing hazard, indicative of aging/degradation). The failure times (in months) are:

   5.2,  7.8,  9.1, 11.3, 12.5, 13.0, 14.2, 15.1, 15.9, 16.7, 17.2, 17.8, 18.4, 18.9, 
  19.3, 19.7, 20.2, 20.6, 21.0, 21.5, 21.9, 22.3, 22.7, 23.1, 23.5, 23.9, 24.3, 24.7, 
  25.1, 25.5, 25.9, 26.3, 26.7, 27.1, 27.5, 27.9, 28.3, 28.7, 29.1, 29.5, 29.9, 30.3, 
  30.7, 31.1, 31.5, 31.9, 32.3, 32.7, 33.1, 33.5, 33.9, 34.3, 34.7, 35.1, 35.5, 35.9, 
  36.3, 36.7, 37.1, 37.5, 37.9, 38.3, 38.7, 39.1, 39.5, 39.9, 40.3, 40.7, 41.1, 41.5,
  41.9, 42.3, 42.7, 43.1, 43.5

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

We test

\(H_0: \beta = 1\) (Exponential model)

versus

\(H_1: \beta \ne 1\) (Weibull model).

Under \(H_0\), the hazard is constant. Under \(H_1\), the hazard may increase or decrease with time.

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 and \(\beta>0\) is the shape parameter. so that \[ \ell(\lambda, \beta) = \sum_{i=1}^{n} \log f(x_i; \lambda, \beta) = n \log \beta - n \beta \log \lambda + (\beta - 1) \sum_{i=1}^{n} \log x_i - \sum_{i=1}^{n} \left( \frac{x_i}{\lambda} \right)^\beta \]

The gradient functions (partial derivatives) are \[ \frac{\partial \ell}{\partial \lambda} = \frac{\beta}{\lambda} \left[ \sum_{i=1}^n \left( \frac{x_i}{\lambda} \right)^\beta - n \right] \] \[ \frac{\partial \ell}{\partial \beta} = \frac{n}{\beta} - n \log \lambda + \sum_{i=1}^{n} \log x_i - \sum_{i=1}^{n} \left( \frac{x_i}{\lambda} \right)^\beta \log \left( \frac{x_i}{\lambda} \right) \]

x <- c(
  5.2, 7.8, 9.1, 11.3, 12.5, 13.0, 14.2, 15.1, 15.9, 16.7, 17.2, 17.8, 18.4, 18.9,
  19.3, 19.7, 20.2, 20.6, 21.0, 21.5, 21.9, 22.3, 22.7, 23.1, 23.5, 23.9, 24.3, 24.7,
  25.1, 25.5, 25.9, 26.3, 26.7, 27.1, 27.5, 27.9, 28.3, 28.7, 29.1, 29.5, 29.9, 30.3,
  30.7, 31.1, 31.5, 31.9, 32.3, 32.7, 33.1, 33.5, 33.9, 34.3, 34.7, 35.1, 35.5, 35.9,
  36.3, 36.7, 37.1, 37.5, 37.9, 38.3, 38.7, 39.1, 39.5, 39.9, 40.3, 40.7, 41.1, 41.5,
  41.9, 42.3, 42.7, 43.1, 43.5
)
n <- length(x)

# Log-likelihood function
weibull_loglik <- function(par, data) {

  lambda <- par[1]
  beta   <- par[2]

  # Parameter constraint
  if (lambda <= 0 || beta <= 0) return(-Inf)

  ll <- length(data) * log(beta) -
        length(data) * beta * log(lambda) +
        (beta - 1) * sum(log(data)) -
        sum((data / lambda)^beta)

  return(ll)
}

# Negative log-likelihood (for minimization)
weibull_nll <- function(par, data) {
  return(-weibull_loglik(par, data))
}

# Gradient function
weibull_grad <- function(par, data) {

  lambda <- par[1]
  beta   <- par[2]
  n <- length(data)

  z <- (data / lambda)^beta

  # derivative wrt lambda
  d_lambda <- - (beta / lambda) * (sum(z) - n)

  # derivative wrt beta
  d_beta <- - (n / beta -
               n * log(lambda) +
               sum(log(data)) -
               sum(z * log(data / lambda)))

  return(c(d_lambda, d_beta))
}

# Run optimization
fit_weibull <- optim(
  par = c(mean(x), 1.5),   # initial values
  fn  = weibull_nll,
  gr  = weibull_grad,
  data = x,
  method = "L-BFGS-B",
  lower = c(1e-8, 1e-8)
)

fit_weibull$par

[1] 31.418211  3.370789

Summary: Using optim(), the maximum likelihood estimates are \(\hat{\beta} \approx 3.37\). Thus the fitted Weibull model is \(\hat{\lambda} \approx 31.42\).

Since \(\hat{\beta} > 1\), the fitted model suggests an increasing hazard rate, which is consistent with mechanical wear and aging.

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

Step 1: Model

When \(\beta =1\) \(f(t;\lambda) = \frac{1}{\lambda} e^{-t/\lambda}\)

Step 2: Log likelihood

\(\ell(\lambda) = -n \log \lambda - \frac{1}{\lambda} \sum_{i=1}^{n} x_i\)

Step 3: Derivative

\(\frac{d\ell}{d\lambda} = -\frac{n}{\lambda} + \frac{\sum x_i}{\lambda^2}\) Set to zero:\(\hat{\lambda} = \bar{x}\)

# Exponential MLE
lambda0_hat <- mean(x)

lambda0_hat

[1] 28.18533

Summary:

The estimate is\(\hat{\lambda}_0 = \bar{x} \approx 28.1853\).

c). Use a) and b) to perform the regular likelihood ratio \(\chi^2\) test for \(\beta = 1\) and report the p-value.

Step 1: Statistic

\(LR = -2[\ell_0 - \ell_1]\)

Step 2:

# Exponential log-likelihood
exp_loglik <- function(lambda, data) {
  -length(data) * log(lambda) - sum(data) / lambda
}

# Compute log-likelihoods
logLik0 <- exp_loglik(lambda0_hat, x)

lambda1_hat <- fit_weibull$par[1]
beta1_hat   <- fit_weibull$par[2]

logLik1 <- weibull_loglik(c(lambda1_hat, beta1_hat), x)

# Likelihood ratio
LR_obs <- -2 * (logLik0 - logLik1)

# p-value (chi-square)
p_chisq <- 1 - pchisq(LR_obs, df = 1)

LR_obs

[1] 100.6144

p_chisq

[1] 0

Summary: Because \(LR \approx 100.61\) \(p < 2.2 \times 10^{-16}\) Therefore, Reject \(H_0\).

d). Use the BLRT algorithm to perform a bootstrap likelihood ratio test and report the bootstrap p-value. Note that you are expected to translate the BLRT algorithm into R code to perform the BLRT. [Hint: The chi-square distribution should not be used in this part of the analysis.]

Step 1: IdeaInstead of: \(LR \sim \chi^2\) we approximate the distribution using bootstrap.

Step 2: Algorithm Fit model under \(H_0\). Then, generate bootstrap samples from exponential and compute LR each time. Last, Compare with observed.

Step 3:

set.seed(123)

# Fit exponential
fit_exp_mle <- function(data) {
  lambda_hat <- mean(data)
  logLik <- -length(data) * log(lambda_hat) - sum(data) / lambda_hat
  return(list(lambda_hat = lambda_hat, logLik = logLik))
}

# Fit Weibull
fit_weibull_mle <- function(data) {
  out <- optim(
    par = c(mean(data), 1.5),
    fn = weibull_nll,
    gr = weibull_grad,
    data = data,
    method = "L-BFGS-B",
    lower = c(1e-8, 1e-8)
  )

  return(list(
    lambda_hat = out$par[1],
    beta_hat   = out$par[2],
    logLik     = -out$value
  ))
}

# Observed statistic
fit0 <- fit_exp_mle(x)
fit1 <- fit_weibull_mle(x)

LR_obs <- -2 * (fit0$logLik - fit1$logLik)

# Bootstrap
B <- 2000
LR_boot <- numeric(B)

for (b in 1:B) {

  # Generate from H0
  x_star <- rexp(n, rate = 1 / fit0$lambda_hat)

  fit0_b <- fit_exp_mle(x_star)
  fit1_b <- fit_weibull_mle(x_star)

  LR_boot[b] <- -2 * (fit0_b$logLik - fit1_b$logLik)
}

# Bootstrap p-value
p_boot <- (1 + sum(LR_boot >= LR_obs)) / (B + 1)

p_boot

[1] 0.0004997501

Summary:

Using \(B=2000\) bootstrap samples, the bootstrap p-value is:\(p_{boot} \approx 0.0005\)This is extremely small, so the BLRT also rejects \(H_0\).Because bootstrap results depend slightly on the random seed and the chosen \(B\), the exact value may vary a little, but it will still be very small.DecisionSince \(p_{boot} \approx 0.0005 < 0.05\), we again reject \(H_0: \beta = 1\).

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

• Whether the two tests generated the same results.

Summary: Both the classical likelihood ratio test and the bootstrap LRT reject \(H_0: \beta=1\), with p-values near zero (\(p < 2.2 \times 10^{-16}\) and \(p_{boot} \approx 0.0005\)). This provides strong evidence that the exponential model does not fit the data. The estimated Weibull shape parameter \(\hat{\beta} \approx 3.37 > 1\) indicates an increasing hazard rate over time, consistent with mechanical wear. Therefore, the Weibull model is preferred, while the exponential model underfits the data.

• Which model is recommended for the data.

Summary: Weibull model is recommended because it captures increasing failure risk and fits the data significantly better.