STA 506 Homework 6: Concepts and Construction of Confidence Intervals

1 Introduction

For this assignment, we will be estimating the scale parameter of the Lindley distribution, \(\theta\), by using confidence intervals. Specifically, we will build a 95% asymptotic confidence interval and a 95% likelihood ratio confidence interval to estimate \(\theta\).

In order to build the confidence intervals, we will be using a data set consisting of customer service call times (in minutes) from a major telecommunications provider in North America.

calls <- c(3.2, 5.8, 7.1, 4.5, 10.3, 6.2, 8.7, 5.1, 12.5, 6.9,
9.4, 5.7, 11.8, 4.9, 9.1, 6.5, 13.2, 7.8, 10.6, 6.1,
8.9, 5.4, 12.1, 7.3, 9.8, 5.9, 11.4, 6.8, 10.9, 7.5,
4.2, 8.3, 6.4, 14.1, 5.6, 9.7, 7.9, 11.1, 6.7, 10.2,
5.3, 8.6, 7.2, 12.9, 6.3, 9.3, 8.1, 13.7, 7.6, 10.8)

2 Asymptotic Confidence Interval

We will first build the asymptotic confidence interval. In order to build an asymptotic confidence interval, we will need to check if we can use the central limit theorem. According to the central limit theorem, the distribution of a sample mean will approach a normal distribution as the sample size increases, regardless of what the initial distribution looked like. In particular, the central limit theorem can be used if the sample size is greater than or equal to 30.

n <- length(calls) # Sample Size
print(paste("Sample Size:", n))

[1] "Sample Size: 50"

As we can see, our data set has 50 observations, so the central limit can apply. So we can approximate \(\theta\) and construct our asymptotic confidence interval using the normal distribution.

To start, we can find a point estimate for \(\theta\) using the maximum likelihood method. The maximum likelhood estimator for the lindley distribution is:

\[ \hat{\theta} = \frac{1 - \bar{x} + \sqrt{\bar{x}^2 + 6\bar{x} + 1}}{2\bar{x}} \]

Using this formula, we obtain a maximum likelihood estimate of 0.2191.

mean.calls <- mean(calls) # Sample Mean

theta.mle <- (1 - mean.calls + sqrt((mean.calls^2) + 6*mean.calls + 1))/(2*mean.calls) # Maximum Likelihood Estimate for Theta
theta.mle

[1] 0.2190994

Next, we need the margin of error. The margin of error consists of two parts: a critical value and the standard error. The critical value can be obtained via the qnorm function.

# Calculating critical value using alpha = 0.05/2 due to symmetry
critvalue.95 <- qnorm(1 - 0.05/2, mean = 0, sd = 1, lower.tail = TRUE) 
critvalue.95

[1] 1.959964

The standard error in this case will be the sqrt of the inverse of the fisher information. The fisher information of the lindley distribution is: \[ I(\theta) = \frac{2}{\theta^2} - \frac{1}{(1+\theta)^2} \]

Using this formula, we obtain a fisher information of 7.4884.

# Formula for Fisher Information of Lindley Distribution
fish.info <- (2/(theta.mle^2)) - (1/(1 + (theta.mle^2)))

fish.info

[1] 40.70853

Finally, multiplying the critical value and the standard error will give us the margin of error.

var.mle <- 1/n*fish.info

se.mle <- sqrt(var.mle)
se.mle

[1] 0.9023141

moe <- critvalue.95*se.mle
moe

[1] 1.768503

ci.lowlim.asymp <- theta.mle - moe # Lower Limit of Confidence Interval
ci.uplim.asymp <- theta.mle + moe # Upper Limit of Confidence Interval
con.int.asymp = c(ci.lowlim.asymp, ci.uplim.asymp) # Combined Confidence Interval
con.int.asymp

[1] -1.549404  1.987602

3 Likelihood Ratio Confidence Interval

Next, we will build the likelihood ratio confidence interval. The likelihood ratio confidence interval are more versatile and more accurate ways of constructing confidence intervals, especially if the sample size is not sufficiently large.

Our likelihood confidence interval is the range of values that satisfies Wilks’ theorem.

Wilks’ Theorem: Assume that \(\theta\) is the true unknown population parameter If sample size \(n\) is large, we have

\[ LR(\theta)=-2\log\Lambda(\theta) = 2\left[ \log L(\hat{\theta}_{MLE}; x) -\log L(\theta; x)\right]\xrightarrow{d} \chi^2_1 \]

For simplicity of calculations, instead of using the standard likelihood function, we will use the log-likelihood function.

\[ L(\theta) = \prod_{i=1}^n f(x_i;\theta) = \prod_{i=1}^n \left[ \frac{\theta^2}{1+\theta} (1 + x_i) e^{-\theta x_i} \right]. \]

Let \(S = \sum_{i=1}^n x_i\), \(\bar{x} = S/n\), and \(C = \sum_{i=1}^n \ln(1 + x_i)\) (constant with respect to \(\theta\)):

\[ \ell(\theta) = \ln L(\theta) = n \ln\left( \frac{\theta^2}{1+\theta} \right) + C - \theta S. \]

To start, we need a point estimate for theta. Since this is a likelihood ratio confidence interval, we will need the maximum likelihood estimator. We already established when constructing the asymptotic interval that our mle is \(\hat{theta}\) = 0.2191.

theta.mle # Our Maximum Likelihood Estimate from Before

[1] 0.2190994

Next, we can obtain our critical value via the qchisq() function.

# Critical Chi-Square Value for Confidence Interval
crit.chisqr <- qchisq(0.95, df = 1, lower.tail = TRUE) 
crit.chisqr

[1] 3.841459

Next we can obtain our maximum log-likelihood estimate:

# Defining Log-Likelihood Function
log.lik <- function(theta, data){
n = length(data)
c = sum(log(1 + data))
s = sum(data)
  
  return(n*log(theta^2/(1+theta)) + c - (theta * s)) 
  # Returns log likelihood estimate
}

log.mle <- log.lik(theta.mle, calls) # function evaulated at MLE
log.mle

[1] -143.3101

We can now utitlize Wilks’ Theorem to find our bounds for our confidence interval:

# Utilizing Wilks' Theorem
Wilks <-  function(theta, data){
  2*(log.mle - log.lik(theta, data))
}

roots <- function(theta, data, theta.mle){
  log.mle - crit.chisqr
}