Assignment 7: Constructing Likelihood Ratio Confidence Interval

Assignment Objectives

• Reinforce the likelihood concepts and MLE.

• Understand the concepts of confidence intervals.

• Master the process of finding likelihood ratio confidence interval of unknown parameter.

Policies of Using AI Tools

Policy on AI Tool Use: You must 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.

One Parameter Lindley Distribution

The Lindley distribution is a continuous probability distribution proposed by D.V. Lindley in 1958. It represents a weighted mixture of exponential and gamma distributions, providing a flexible single-parameter model for lifetime data.

\[ f(x;\theta) = \frac{\theta^2}{1+\theta}(1+x)e^{-\theta x}, \quad x > 0, \quad \theta > 0 \]

where \(x\) = random variable (e.g., time, size, amount) and \(\theta\) = shape parameter controlling the distribution.

Given an independent random sample \(X_1, X_2, \dots, X_n\):

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

After some algebra, we obtain the closed form of the MLE of \(\theta\) in the following

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

As good exercise, we can derive the following Fisher information of \(\theta\):

\[ \boxed{I(\theta) = \frac{2}{\theta^2} - \frac{1}{(1+\theta)^2}} \]

This assignment focuses on constructing various confidence intervals of the shape parameter \(\theta\) in the Lindley distribution.

Question: Customer Service Times (minutes)

The customer service call duration data set originates from a major telecommunications provider in North America, operating in a highly competitive market where:

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

Assuming the data follow a one-parameter Lindley distribution, construct a \(95\%\) confidence interval for the parameter \(\theta\) using the provided data and the specified methods. For each of the following questions, first describe your reasoning process for the analysis, then write code to perform the actual analysis. Finally, summarize the results to conclude the question.

1. Construct a 95% asymptotic confidence interval based on the asymptotic sampling distribution of the maximum likelihood estimator (MLE) of \(\theta\).

Question 1 Part A

First we will mathematically determine the formula for calculating the 95% asymptotic confidence interval based on our formulas for \(\hat{\theta}\) and \(I(\theta)\). Given the large population size we and finite variance we can reasonably assume:

\[ \hat{\theta} \xrightarrow{d} N(\theta, \widehat{Var}(\hat{\theta})) \]

We can then standardize our formula:

\[ \frac{\hat{\theta}-\theta}{\sqrt{\widehat{Var}(\hat{\theta})}}\xrightarrow{d} N(0,1) \]

Finally, we can construct the 95% confidence interval:

\[ \hat{\theta} \pm Z_{1-0.05/2}{\sqrt{\widehat{Var}(\hat{\theta})}} \]

Since this is normally distributed, we can use the derived Fisher information to find the variance:

\[ I(\theta)=\frac{1}{n(Var(\theta))} \Rightarrow Var(\theta)=\frac{1}{n(I(\theta))} \]

Where \(I(\theta) = \frac{2}{\theta^2} - \frac{1}{(1+\theta)^2}\).

Finally, we can use the closed form MLE above to find \(\hat{\theta}\):

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

We can copy these formulas into our code to perform the analysis on our dataset:

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

n <- length(tele) #used to get n
m <- mean(tele) #used to find mean
theta.hat <- (1-m+sqrt(m**2+6*m+1))/(2*m) #formula for theta-hat
i <- 2/(theta.hat**2)-1/((1+theta.hat)**2) #formula for i(theta)
se <- sqrt(1/(n*i)) #formula for standard deviation

ci.low <- theta.hat - qnorm(1-0.025)*se #calculating lower limit
ci.up <- theta.hat + qnorm(1-0.025)*se #calculating upper limit
ci <- c(ci.low, ci.up)

cat("95% Asymptotic Confidence Interval for Theta-Hat: \n",
    "Asymptotic CI: [", ci, "] \n") #stating interval

95% Asymptotic Confidence Interval for Theta-Hat: 
 Asymptotic CI: [ 0.1758057 0.2623931 ] 

In conclusion, we are 95% confident that the true value of \(\theta\) for the population is between 0.1758057 and 0.2623931.

2. Construct a 95% likelihood ratio confidence interval for \(\theta\).

Question 1 Part B

From class, we know that the likelihood ratio statistic for this setting would be:

\[ 2(\ell(\hat{\theta})-\ell(\theta)) \le \chi_{1,0.95}^2 \] Which means that the interval endpoints \(\theta_0\) and \(\theta_1\) of rhe \(100(1-\alpha)%\) CI should satisfy:

\[ \ell(\theta)\ge-\ell(\hat{\theta})-\chi_{1,0.95}^2/2 \] Furthermore, we know that the formula for \(\ell(\theta)\) is:

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

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

And the formula for \(\hat{\theta}\) is:

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

We can then use uniroot() solve this for our confidence interval.

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

n <- length(tele) #used to get n
sumx <- sum(tele) #used to get sum of obs
m <- mean(tele) #used to get mean
theta.hat <- (1-m +sqrt(m^2+6*m+1))/(2*m) #formula for theta.hat
loglik <- function(theta) { #l(theta) formula
  n*log(theta**2/(1+theta))+sum(log(1+tele))-theta*sumx
}
lr <- function(theta) { #likelihood ration function
  2*(loglik(theta.hat)-loglik(theta))-qchisq(0.95, 1)
}

ci.low2 <- uniroot(lr,lower=1e-6, upper=theta.hat)$root #Used to get bounds
ci.up2 <- uniroot(lr,lower=theta.hat, upper=5)$root

cat("95% Likelihood Ratio CI: [", ci.low2, ",", ci.up2, "]\n") #states bounds

95% Likelihood Ratio CI: [ 0.1786411 , 0.2653396 ]

In conclusion, we are 95% confident that the true value of \(\theta\) for the population is between 0.1786411 and 0.2653396.

3. Assuming the two confidence intervals above are valid, compare them in terms of performance and make a recommendation. Justify your recommendation.

Answer to Question 1 Part C

Assuming that both of the above confidence intervals are valid, in that for both of them we can be at least 95% confident the true value of \(\theta\) is between their lower and upper limits, I would recommend using the asymptotic approach given that it resulted in a slightly smaller interval of 0.0865874 in length compared to likelihood ratio confidence interval which was 0.0866985 in length. This is a small difference, but if I was making a recommendation between the two I would recommend the one that had a slightly smaller interval.