STA 506 Assignment 10: Wald, Score, and Likelihood Ratio Tests

1 Introduction

For this assignment, we will focus on performing asymptotic \(\chi^2\) hypothesis tests to determine whether or not a sample came from one distribution or another. Specifically, we will do a likelihood ratio and Wald \(\chi^2\) to determine if a sample came from a weibull distribution or an exponential distribution.

The sample we will use for our tests is a sample consisting of failure times (in hours) for bearings used in motor assemblies at a distribution company’s facility.

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

2 Part A:

For this part, we will find the maximum likelihood estimators of the Weibull distribution, treating this sample as if it comes from a Weibull distribution.

First, we’ll start with the log-likelihood function for the weibull distribution: $$ \[\begin{align} \ell(k, \lambda; \mathbf{t}) &= \ln L(k, \lambda; \mathbf{t}) \\ &= \sum_{i=1}^{n} \ln \left[ \frac{k}{\lambda} \left( \frac{t_i}{\lambda} \right)^{k-1} \exp\left( -\left( \frac{t_i}{\lambda} \right)^k \right) \right] \\ &= \sum_{i=1}^{n} \left[ \ln k - \ln \lambda + (k-1)(\ln t_i - \ln \lambda) - \left( \frac{t_i}{\lambda} \right)^k \right] \\ &= \sum_{i=1}^{n} \left[ \ln k - k \ln \lambda + (k-1) \ln t_i - \left( \frac{t_i}{\lambda} \right)^k \right] \end{align}\]

$$

From there, we can take the partial derivatives (one with respect to each variable), set them both equal to 0, and then solve for the two parameters and get their maximum likelihood estimates. If done right, the maximum likelihood estimates are: \[ \begin{align} \hat{\lambda} &= \left( \frac{1}{n} \sum_{i=1}^{n} x_i^{\hat{\beta}} \right)^{1/\hat{\beta}} \quad \\ \frac{1}{\hat{\beta}} &= \frac{\sum_{i=1}^{n} (t_i^{\hat{\beta}} \ln x_i)}{\sum_{i=1}^{n} x_i^{\hat{\beta}}} - \frac{1}{n} \sum_{i=1}^{n} \ln x_i \quad \end{align} \]

Where \(\lambda\) is the scale parameter and \(\beta\) is the shape parameter.

To obtain these parameters, we will use R’s built in optim() function, which will solve a system of equations numerically. The following code will be used to obtain the parameter estimates:

## Manual MLE implementation for Weibull distribution

# Create Weibull log-likelihood function
weibull.loglik <- function(params, data) {
  shape <- params[1]  # passing the parameters
  scale <- params[2]

n <- length(data)   # sample size
  
# Log-likelihood for Weibull distribution
loglik <- n * log(shape) - n * shape * log(scale) + 
            (shape - 1) * sum(log(data)) - 
            sum((data / scale)^shape)
  
  return(loglik)    # Return
}


# Score equation

weibull.score <- function(params, data) {
  shape <- params[1]
  scale <- params[2]
  n <- length(data)
  
# Gradient for shape parameter
grad_shape <- n/shape - n * log(scale) + 
                sum(log(data)) - 
                sum((data/scale)^shape * log(data/scale))
  
# Gradient for scale parameter
grad_scale <- -(n * shape)/scale + 
                (shape/scale) * sum((data/scale)^shape)
  
  return(c(grad_shape, grad_scale))
}

# Need to provide initial values for parameters
initial.params <- c(shape = 1, scale = 5)  # Reasonable starting values


# Using optim with Nelder-Mead method
mle.result.weibull <- optim(
  par = initial.params,
  fn = weibull.loglik,
  gr = weibull.score,
  data = times,
  method = "L-BFGS-B",
  hessian = TRUE,
  control = list(trace = FALSE,
                 fnscale = -1,
                 maxit = 500,
                 abstol = 1e-8)
)
##
mle.result.weibull$par

    shape     scale 
 2.205999 99.020285 

mle.result.weibull$hessian

            shape       scale
shape -17.5910623  0.20545481
scale   0.2054548 -0.02481605

The estimate for the shape parameter is 2.206, while the estimate for the scale parameter is 99.0203.

To verify, we can use the fitdistr() function from the MASS package

mle.weibull.fit <- fitdistr(times, "weibull")
mle.weibull.fit

     shape        scale   
   2.2064616   99.0593897 
 ( 0.2508946) ( 6.6816324)

They’re just about the same, so we can say that our original estimates are correct.

3 Part B

For this part, we will find the maximum likelihood estimator for the exponential distribution, treating this sample as if it came from an exponential distribution.

In the case of the exponential distribution, we are only estimating the scale parameter because an exponential random variable is just a Weibull random variable with the shape parameter, \(\beta\), set to 1.

Depending on how we write the density function for the exponential, the resulting maximum likelihood estimate for the scale parameter is: \[ \hat{\lambda} = \frac{1}{n} \sum x_i = \bar{x} \implies \frac{1}{\hat{\lambda}} = \frac{n}{\sum x_i} = \frac{1}{\bar{x}} \]

Regardless of how the estimator is written, the process for solving for \(\lambda\) is the same. In this case, since we only have one parameter to solve for, and it has a simple form, we can use simpler methods to solve for it. Since we already determined that \(\lambda\) is just the sample mean, we will use that estimate.

First, since we did it manually in Part A, we can do it manually in this part as well.

## Manual MLE implementation for Exponential distribution

# Create Exponential Log-Likelihood Function
exponential.loglik <- function(params, data) {
  scale <- params[1]

n <- length(data)   # sample size
  
# Log-likelihood for Exponential Distribution
loglik <- -n*log(scale) - sum(data)/scale
  
  return(loglik)    # Return 
}


# Score equation

exponential.score <- function(params, data) {
  scale <- params[1]
  n <- length(data)

  
# Gradient for scale parameter
grad_scale <- -(n/scale) + sum(data)/scale^2

  return(grad_scale)
}

# Need to provide initial values for parameters
  # Reasonable starting values


# Using optim with Nelder-Mead method
mle.result.exponential <- optim(
  par = 5, # Just a generic number for an estimate
  fn = exponential.loglik,
  gr = exponential.score,
  data = times,
  method = "L-BFGS-B",
  hessian = TRUE,
  control = list(trace = FALSE,
                 fnscale = -1,
                 maxit = 500,
                 abstol = 1e-8)
)
##
mle.result.exponential$par

[1] 87.61

Now let’s plug the sample mean for the mle:

## Using Sample Mean as MLE For Exponential 

mle.exponential.scale <- mean(times)
mle.exponential.scale

[1] 87.61

The mean of our sample is 87.61, so the estimate for our scale parameter is 87.61.

To verify, we can use fitdistr() from the MASS package:

## Verifying MLE For Exponential via fitdistr()

exponential.fit <- fitdistr(times, "exponential")
exponential.fit

      rate    
  0.011414222 
 (0.001614215)

As we can see, R gave us the estimate for the rate parameter. But we need the estimate for the scale parameter. The good news is that we can get the estimate for the scale parameter, as it is just the inverse of the rate parameter. In other words, scale = 1/rate.

mle.fit <- 1/exponential.fit$estimate
mle.fit

 rate 
87.61 

4 Part C:

For this part, we will conduct the first of two hypothesis tests to determine whether or not this data comes from a Weibull distribution or an exponential distribution. Specifically, we will conduct a likelihood ratio \(\chi^2\) test. Our null hypothesis, \(H_0\), is that this data comes from an exponential distribution. In other words, our null hypothesis is that \(\beta = 1\). Our alternative hypothesis, \(H_a\), is that it our data comes from a Weibull distribution. In other words, our alternative hypothesis is that \(\beta \neq 1\). To put into symbols: \[ H_0: \beta = 1 (Exponential) \\ H_a: \beta \neq 1 (Weibull) \]

This test will be done at significance level \(\alpha = 0.05\).

## Likelihood Ratio Test using both Log-Likelihood Functions

# Weibull Log-Likelihood Function
loglike.weibull <- function(data, lambda, beta){
  n <- length(data)
  n*log(beta)-n*beta*log(lambda)+(beta-1)*sum(log(data))-sum((data/lambda)^beta)
}

# Exponential Log-Likelihood Function
loglike.exponential <- function(data, lambda){ 
  n <- length(data)
  -n*log(lambda)-sum(data)/lambda
}

## Evaluate Log Likelihood @ MLE
scale.est.weibull <- mle.result.weibull$par[[2]]
shape.est.weibull <- mle.result.weibull$par[[1]]
scale.est.exponential <- mle.fit

# Exponential vs Weibull Hypotheses
LogLike.Alt <- loglike.weibull(times, scale.est.weibull, shape.est.weibull)
LogLike.Null <- loglike.exponential(times, scale.est.exponential)

ratio <- 2*(LogLike.Alt - LogLike.Null) #Computes likelihood ratio
ratio

    rate 
34.44995 

## Getting p-value for ratio test
pvalue.ratio <- 1 - pchisq(ratio, df = 1) #Gets p-value
pvalue.ratio

        rate 
4.373531e-09 

Our \(\chi^2\) statistic is 34.44995, and the p-value is 4.37 x \(10^-9\), which is obviously less than 0.05, so it’s safe to say that we can reject the null hypothesis and conclude that this sample came from a Weibull distribution.

5 Part D

For this part, we will conduct the second of our two hypothesis tests to determine whether or not our data came from an exponential or Weibull distribution. Specifically, we will conduct a Wald \(\chi^2\). The null hypothesis, the alternative hypothesis, and the level of significance are all the same as the ones from the likelihood-ratio test from Part C.

First, we have to get the Fisher Information matrix to get the standard error required for the test.

## Obtaining Hessian and Fisher Matrices

result <- optim(
          par = c(1,1), 
          #need to provide initial values to start the iteration.
                                    # Choosing appropriate initial values is critical.
          fn = weibull.loglik,     # Caution: need negative log-likelihood
          gr = weibull.score, 
          data = times, # also nee negative score
          method = "BFGS",          # calling BFGS algorithm
          hessian = TRUE,           # return Hessian matrix
          control = list(           # some controls in numerical iteration
                fnscale = -1,       # turn the minimization to maximization problem
                maxit = 1000,       # stopping rule 1: cap number of iterations
                reltol = 1e-8,      # stopping rule: precision control
                REPORT = 10)         # The algorithm will print progress updates 
                           # every 10 iterations
)

Hess <- result$hessian
Hess

           [,1]        [,2]
[1,] -17.591017  0.20545403
[2,]   0.205454 -0.02481605

Fisher <- -Hess
Fisher

          [,1]        [,2]
[1,] 17.591017 -0.20545403
[2,] -0.205454  0.02481605

Covar <- Fisher^-1
Covar

            [,1]      [,2]
[1,]  0.05684719 -4.867269
[2,] -4.86726882 40.296494

Next, we

## Obtaining Standard Error for the Wald Test

var.shape <- Covar[2, 2] # Gets variance of beta
se.shape <- sqrt(var.shape) # Gets standard deviation of beta

shape.est <- result$par[1] #Gets MLE of beta
shape.est

[1] 2.206002

Wald <- (shape.est - 1)^2/var.shape #Gets Wald test statistic
Wald

[1] 0.03609346

pvalue.wald <- 1 - pchisq(Wald, df = 1)
pvalue.ratio

        rate 
4.373531e-09 

Our Wald chi-square statistic is 0.0361, and the p-value is 4.37 x 10^-9, which is way less than 0.05, which means we can reject our null hypothesis and conclude that our data came from a weibull distribution.

6 Part E

With both of our tests, we rejected the null hypothesis, meaning we concluded that there is evidence to suggest that our data came from a weibull distribution with both tests. Since that is the case, the weibull distribution is recommended for this sample.