Assignment 11: Detecting Overfitting and Overfitting Issues

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

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

Answer to Question 1 Part A

We are provided with the probability density function (PDF) of the Weibull distribution:

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

We can use this to find the likelihood function:

\[ L(\lambda, \beta)=\prod_{i=1}^nf(t_i ; \lambda, \beta)=\prod_{i=1}^n \frac{\beta}{\lambda} \left( \frac{t_i}{\lambda} \right)^{\beta-1} \exp\left[ -\left( \frac{t_i}{\lambda} \right)^\beta \right] \]

Therefore the log-likelihood function is:

$$ \[\begin{aligned} \ell(\lambda, \beta) &= \log(L(\lambda, \beta)) \\ &= \sum_{i=1}^n\left[\log\beta-\log\lambda+(\beta-1)\log(t_i/\lambda)-(t_i/\lambda)^\beta\right] \\ &= n\log\beta-n\log\lambda+(\beta-1)\sum_{i=1}^n\log t_i-\sum_{i=1}^n(t_i/\lambda)^\beta \end{aligned}\]

$$

Using the log-likelihood function we can take the partial derivatives of \(\lambda\) and \(\beta\) to find the gradient functions:

$$ =-+_{i=1}^n(t_i/)\

=-n+_{i=1}^{nt_i-{i=1}^nt_i-{i=1}}n(t_i/)^(t_i/) $$

We can then use our log-likelihood function and gradient functions to find \(\hat{\lambda}\) and \(\hat{\beta}\):

bearings <- 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)

loglik_weibull <- function(par, x){ #Used to impute log-likelihood
  lambda <- par[1]
  beta <- par[2]
  if (lambda <= 0 || beta <= 0) return (Inf)
  n <- length(x)
  ll <- n*log(beta)-n*beta*log(lambda)+(beta-1)*sum(log(x))-sum((x/lambda)^beta)
  
  return(-ll)
}

grad_weibull <- function(par, x){ #Used to impute gradient functions
  lambda <- par[1]
  beta <- par[2]
  if (lambda <= 0 || beta <= 0) return (c(Inf, Inf))
  n <- length(x)
  term <- (x/lambda)^beta
  d_lambda <- -n*beta/lambda+(beta/lambda)*sum(term)
  d_beta <- n/beta-n*log(lambda)+sum(log(x))-sum(term*log(x/lambda))
  
  return(-c(d_lambda, d_beta))
}

result <- optim( #Used to compute MLEs
  par=c(1,1),
  fn=loglik_weibull,
  gr=grad_weibull,
  x=bearings,
  method="L-BFGS-B",
  lower=c(1e-6, 1e-6)
)
result

$par
[1] 99.020320  2.206001

$value
[1] 256.4198

$counts
function gradient 
      26       26 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

lambda_hat <- result$par[1] #Gets MLEs from optim function
beta_hat <- result$par[2]

lambda_hat

[1] 99.02032

beta_hat

[1] 2.206001

The MLE of the Weibull parameters are \(\hat{\lambda}=99.02032\) and \(\hat{\beta}=2.206001\).

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

Answer to Question 1 Part B

We are provided with the probability density function (PDF) of the Weibull distribution when \(\beta=1\):

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

We can use this to find the likelihood function:

\[ L(\lambda)=\prod_{i=1}^nf(t; \lambda)=\prod_{i=1}^n\frac{1}{\lambda}e^{-t_i/\lambda} \]

Then we can find the log-likelihood function:

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

We can then find the gradient function and solve for \(\hat{\lambda}\):

$$ \[\begin{aligned} &\frac{\partial\ell}{\partial\lambda}=-\frac{n}{\lambda}+\frac{1}{\lambda^2}\sum_{i=1}^n t_i \\ &\Rightarrow -\frac{n}{\lambda}+\frac{1}{\lambda^2}\sum_{i=1}^n t_i=0 \\ &\Rightarrow -n\lambda+\sum_{i=1}^n t_i =0 \\ &\Rightarrow \lambda=\frac{1}{n}\sum_{i=1}^n t_i \end{aligned}\]

$$

Therefore,

\[ \hat{\lambda}=\frac{1}{n}\sum_{i=1}^n t_i \]

We can use this to find \(\hat{\lambda}\):

lambda_hat2 = mean(bearings) #Gets lambda MLE

lambda_hat2

[1] 87.61

The MLE of the Weibull parameter when \(\beta=1\) 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).]

Answer to Question 1 Part C

ll_weibull <- function(x, lambda, beta){ #Computes Weibull log-likelihood
  n <- length(x)
  n*log(beta)-n*beta*log(lambda)+(beta-1)*sum(log(x))-sum((x/lambda)^beta)
}

ll_exp <- function(x, lambda){ #Computes exponential log-likelihood
  n <- length(x)
  -n*log(lambda)-sum(x)/lambda
}

lambda_hat_weibull <- 99.020320 #Plugs in MLEs
beta_hat_weibull <- 2.206001
lambda_hat_exp <- 87.61

ll1 <- ll_weibull(bearings, lambda_hat_weibull, beta_hat_weibull) #Gets maximized log-likelihoods
ll0 <- ll_exp(bearings, lambda_hat_exp)

LR <- 2*(ll1 - ll0) #Computes likelihood ratio
LR

[1] 34.44995

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

[1] 4.373531e-09

The p-value is 4.373531e-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.]

Answer to Question 1 Part D

result2 <- optim( #Used to get hessian
  par=c(1,1),
  fn=loglik_weibull,
  gr=grad_weibull,
  x=bearings,
  method="L-BFGS-B",
  lower=c(1e-6, 1e-6),
  hessian=TRUE
)
H <- result2$hessian
cov_mat <- solve(H) #Used to get covariance matrix
cov_mat

           [,1]       [,2]
[1,] 44.6100636 0.52102376
[2,]  0.5210238 0.06293246

var_beta <- cov_mat[2,2] #Gets variance of beta
se_beta <- sqrt(var_beta) #Gets standard deviation of beta

beta_hat <- result2$par[2] #Gets MLE of beta
W <- (beta_hat-1)^2 /var_beta #Gets Wald test statistic
W

[1] 23.11111

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

[1] 1.529047e-06

The p-value is 1.529047e-06.

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.

Answer to Question 1 Part E

Both tests indicated the \(\beta \ne 1\) (at the 0.05 significance level) given that \(\text{p-value}= 4.373531 \times 10^{-09} < 0.05\) and \(\text{p-value}= 1.529047 \times 10^{-06} < 0.05\). Therefore, the Weibull model would be recommended for this data given that we have evidence that \(\beta \ne 1\).

x <- seq(min(bearings), max(bearings), length.out=300) #gets x values used later
weib_dens <- (beta_hat_weibull / lambda_hat_weibull) * 
  (x / lambda_hat_weibull)^(beta_hat_weibull - 1) *
  exp(-(x / lambda_hat_weibull)^beta_hat_weibull) #Weibull pdf
exp_dens <- (1 / lambda_hat_exp) * exp(-x / lambda_hat_exp) #Exponential pdf

plot(x, weib_dens, type = "l", lwd = 2, col = "blue", xlab = "", ylab = "Density", main = "Weibull vs. Exponential Densities")
lines(x, exp_dens, lwd = 2, col = "red", lty = 2)
legend("topright",
       legend = c("Weibull (MLE)", "Exponential (MLE)"),
       col = c("blue", "red"),
       lwd = 2,
       lty = c(1, 2)) #Plots densities