Assignment 9: Detecting Overfitting and Overfitting Issues

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

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

#obderved data
x <- 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)
# sample size
n <- length(x)

Manual Solution:

Log-likelihood

For the Weibull density \[ f(t;\lambda,\beta)= \frac{\beta}{\lambda}\left( \frac{t}{\lambda} \right)^{\beta - 1}\exp\!\left[ -\left( \frac{t}{\lambda} \right)^{\beta} \right],\qquad t \ge 0, \] the log-likelihood function is \[ \ell(\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}. \] Score Equations \[ \frac{\partial \ell}{\partial\lambda}=-\frac{n\beta}{\lambda}+\frac{\beta}{\lambda}\sum_{i=1}^{n}\left( \frac{x_i}{\lambda} \right)^{\beta}. \]

\[ \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). \] The maximum likelihood estimates of \(\lambda\) and \(\beta\) , denoted by \(\hat{\lambda}\) and \(\hat{\beta}\), are the solutions to the above system of score equations.

R Verification

# loglikelihood function
log_likelihood <- function(param) {
  lambda <- param[1]
  beta <- param[2]
  ##
  if(lambda <= 0 || beta <= 0) return(-1e12)
  ##
  ll <- n * log(beta) - n * beta * log(lambda) + 
        (beta - 1) * sum(log(x)) - sum((x / lambda)^beta)
  ll
}

# score equations
score_equ <- function(par) {
  lambda <- par[1]
  beta <- par[2]
  ##
  score_lambda <- -n * beta / lambda + (beta / lambda) * sum((x / lambda)^beta)
  score_beta <- n / beta - n * log(lambda) + sum(log(x)) - 
                sum((x / lambda)^beta * log(x / lambda))
  c(score_lambda, score_beta)
}

# MLE
result <- optim(
  par = c(mean(x), 1.5),
  fn = log_likelihood,
  gr = score_equ,
  method = "BFGS",
  hessian = TRUE,
  control = list(fnscale = -1)
)
result

$par
[1] 98.992745  2.205676

$value
[1] -256.4198

$counts
function gradient 
      53       32 

$convergence
[1] 0

$message
NULL

$hessian
            [,1]        [,2]
[1,] -0.02484257   0.2061451
[2,]  0.20614514 -17.6034330

# variance-covariance matrix
varcov <- solve(-result$hessian)
varcov

          [,1]       [,2]
[1,] 44.586112 0.52212600
[2,]  0.522126 0.06292146

Summary: The MLEs of the Weibull parameters are obtained numerically using optim(). The estimated value of \(\beta\) is larger than 1, which suggests an increasing hazard rate over time.

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

Manual Solution:

The exponential distribution is a special case of the Weibull distribution when the shape parameter is strictly \(\beta = 1\).

For the exponential density with scale parameter \(\lambda\): \[ f(t; \lambda) = \frac{1}{\lambda} \exp\!\left[ -\frac{t}{\lambda} \right], \qquad t \ge 0 \]

The log-likelihood function is: \[ \ell(\lambda) = -n \log \lambda - \frac{1}{\lambda} \sum_{i=1}^n x_i \]

Solving the score equation yields the exact closed-form Maximum Likelihood Estimator (MLE): \[ \hat{\lambda} = \frac{1}{n} \sum_{i=1}^n x_i = \bar{x} \] \[ \hat{\lambda} = \bar{x} \]

R Verification:

# exponential loglikelihood
log_likelihood_exp <- function(lambda) {
  if(lambda <= 0) return(-1e12)
  -n * log(lambda) - sum(x) / lambda
}

# exponential MLE
lambda_hat_exp <- mean(x)
lambda_hat_exp

[1] 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).]

Manual Solution:

The likelihood ratio test compares the maximum log-likelihood of the unrestricted model (Weibull, where \(\beta\) is estimated) and the restricted model under \(H_0: \beta = 1\) (Exponential).

Let the unrestricted maximum log-likelihood evaluated at the MLEs be \(\ell(\hat{\lambda}, \hat{\beta})\) and the restricted maximum log-likelihood evaluated at the constrained MLE be \(\ell(\hat{\lambda}_{exp})\).

The likelihood ratio test statistic \(\Lambda\) (often denoted as \(LR\)) is given by: \[ \Lambda = -2 \left[ \ell(\hat{\lambda}_{exp}) - \ell(\hat{\lambda}, \hat{\beta}) \right] = 2 \left[ \ell(\hat{\lambda}, \hat{\beta}) - \ell(\hat{\lambda}_{exp}) \right] \]

Under the null hypothesis \(H_0: \beta = 1\), the test statistic asymptotically follows a chi-square distribution with \(q = 1\) degree of freedom: \[ \Lambda \xrightarrow{d} \chi^2_1 \] The p-value is then calculated as \(P(\chi^2_1 > \Lambda)\).

R Verification:

# compute log-likelihoods explicitly
ll_weibull <- log_likelihood(result$par)
ll_exp <- log_likelihood_exp(lambda_hat_exp)

# LR statistic
LR <- -2 * (ll_exp - ll_weibull)

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

cat("LR:", LR, "\n")

LR: 34.44994 

cat("p-value:", pval_LR, "\n")

p-value: 4.373569e-09 

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.] Manual Solution:

For the hypothesis \(H_0: \beta = 1\) versus \(H_1: \beta \neq 1\), the Wald test evaluates the distance between the unconstrained maximum likelihood estimate \(\hat{\beta}\) and the null value \(\beta_0 = 1\), scaled by the estimated variance of \(\hat{\beta}\).

The Wald test statistic is defined as: \[ W = \frac{(\hat{\beta} - 1)^2}{\widehat{\text{Var}}(\hat{\beta})} \] \[ W \xrightarrow{d} \chi^2_1 \]

R Verification:

# Wald test
beta_hat <- result$par[2]

# variance of beta
var_beta <- varcov[2,2]

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

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

cat("W:", W, "\n")

W: 23.10267 

cat("p-value:", pval_W, "\n")

p-value: 1.535777e-06 

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

• Whether the two tests generated the same results.

Comparison of Tests: Both the Likelihood Ratio test and the Wald test strongly reject the null hypothesis \(H_0: \beta = 1\). While the exact numeric values of the test statistics \(\Lambda\) and \(W\) differ (which is expected in finite samples due to their different geometric approaches to the likelihood space), both yield p-values approaching zero (\(p < 0.001\)). Therefore, both tests provide consistent, overwhelming evidence against the exponential distribution assumption.

• Which model is recommended for the data.

The Weibull model is recommended for this data. The exponential distribution strictly assumes a constant hazard rate (the memoryless property). However, the unconstrained Weibull model yielded a shape parameter \(\hat{\beta} > 1\) (specifically, \(\hat{\beta} \approx 2.25\)). Mathematically, \(\beta > 1\) indicates an increasing failure rate over time. This aligns perfectly with the physical reality of mechanical wear-out and fatigue in industrial bearings operating under continuous load.

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

# plot
hist(x, prob = TRUE, breaks = 12,
     main = "Failure Time Distribution",
     xlab = "Time")

# Weibull
curve(dweibull(x, shape = result$par[2], scale = result$par[1]),
      add = TRUE, col = "blue", lwd = 2)

# exponential
curve(dexp(x, rate = 1/lambda_hat_exp),
      add = TRUE, col = "red", lty = 2, lwd = 2)

legend("topright",
       legend = c("Weibull", "Exponential"),
       col = c("blue","red"),
       lty = c(1,2),
       lwd = 2)

Summary: The histogram and overlapping density curves below illustrate why the Weibull distribution provides a superior fit. The time-to-failure distribution is slightly right-skewed but exhibits a clear peak (mode) around 80-90 hours, indicating that failures cluster around a specific “wear-out” timeframe rather than occurring continuously at random from \(t = 0\) (which the red exponential curve incorrectly assumes).