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\).

• Describe reasoning process for the analysis

We identify the appropriate confidence interval to be approximated with a normal distribution with an assumed large sample size through the asymptotic confidence interval.

Since the variance for the Lindley distribution cannot be directly assumed/found hence being a combined distribution type, we use the Fisher Information to learn more about the information the variance can hold for the Asymptotic Confidence Interval.The use the Fisher information for calculations for the variance and standard error to calculate the confidence intervals.

Mathematical Algorithm

Sort and mean the call times dataset
find the length of the dataset
Calulate the predicted theta parameter for the dataset
Calculate the Fisher Information
Calculate the variance treating the fisher information as the theta of interset for the standardization equation assuming normal distrubution
Construct the confidence interval assuming normal distrubution 

calltimes= sort(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(calltimes)

mean_call = mean(calltimes)

theta_predict = ((1 - mean_call) + sqrt((mean_call^2 + 6*mean_call + 1))) / (2 * mean_call) #Closed form of the MLE for theta

#First we find the asymptotic distribution prediction parameter

#Goal construct CI for theta --> Asymptotic Confidence Intervals

##Calculating the fisher information for MLE theta

##---
###Building the denominator of the STANDARDIZE (theta_hat - theta)/  sqrt(hat.var(theta))

fisherInfo= (2/theta_predict^2)-(1/ (1+ theta_predict)^2) # a measurement to help measure variance 



#Calculate the variance of the estimator
theta_variance =1/ (n*fisherInfo)

#Calculate standard error (sqrt(varience))

se_theta = sqrt(theta_variance)

#Calculating the Asymptotic CIs
#Construct CI
#z crit is set to 0.975 bc alpha =0.05 --> (100%-(5%/2) = 97.5% =0.975

z = qnorm(0.975)
lower_ci= theta_predict - z * se_theta
upper_ci = theta_predict + z * se_theta

asym_bounds=c(lower_ci, upper_ci)
width_asym = upper_ci -lower_ci

cat("Asymptotic 95% Confidence Interval: [", asym_bounds, "]")

Asymptotic 95% Confidence Interval: [ 0.1758057 0.2623931 ]

#---------------------------------try
#Find asymptotic distribution:
#qnorm()


#Standardize

#(theta_predict-calltimes)/sqrt(var(calltimes))


#Construct CI
#z crit is set to 0.975 bc alpha =0.05 --> (100%-(5%/2) = 97.5% =0.975
#z = qnorm(0.975)

#theta_predict +

Summarize results for question 1a \ The MLE parameter \(\theta\) = 0.2190994has a asymptotic 95%CI [0.1758057, 0.2623931]. The 95%CI bounds include the predicted parameter within its bounds nor do the bounds contain 0 suggesting some strength to the value of ranges the MLE parameter can take on. The difference between the confidence interval ( 0.2623931 - 0.1758057) has an width of 0.0865874, which by the naked eye can appear to be a small number suggesting stronger relevance with statistical significance. Using the asymptotic distribution the telemarketer call time is estimated to be 0.2190994 and we are 95% confident the true call time is between [0.1758057, 0.2623931].

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

• Describe reasoning process for the analysis

We use likelihood ratio confidence intervals to construct confidence intervals for parameters in parametric models. Likelihood ratio confidence intervals present strong accurate interval suited for small and moderate sample sizes and are even more accurate than asymptotic normality confidence intervals.

Mathematical Algorithm

Sort and mean of the call times dataset
Find the number of (length) of values in the call times dataset
Calculate the predicted parameter value for the MLE
-In most cases find the derivative in respect to the parameter (not needed with example given)
Find the maximum log likelihood
Use the maximum log likelihood to find the likelihood ratio statistic
Use the chi square distrubution as quantiled to estimate the interval endpoints to define the confidence intervals for theta
Define the 95%CI for the Lindley (Call times) Distrubution 

calltimes= sort(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(calltimes)

mean_call = mean(calltimes)

#Find the MLE
theta_predict = ((1 - mean_call) + sqrt((mean_call^2 + 6*mean_call + 1))) / (2 * mean_call) 

#we dont need to take the derivative in respect to theta & solve for theta
#to find the log likelihood function directly


#Compute the l_max from the log of the parameter

#The maximum log-likelihood:
l_max = n * log(theta_predict^2/(1+theta_predict)) + sum(log(1+calltimes)) -(theta_predict*sum(calltimes))

#Since we have the log likelihood we can find the log likelihood ratio

###Theta_like_ratio = 2 * (l_max - log(theta_predict))  ##REPLACE THIS WITH CODE

#Find the expected crit value for chi square distribution 
crit = qchisq(0.95, 1) / 2 

thres_height= l_max- crit





####################################################################
log_lik_lindley= function(theta, data)
{
  n = length(data)
  S = sum(data)
  C = sum(log(1+data)) #constant term
  return( 
    n * log(theta^2/(1+theta)) +sum(log(1 + data)) - (theta * sum(data))
    
    )

  
}


#LR_statistic = 2 * (l_max - log_lik_lindley(theta, calltimes))

#Find bounds and check gaps between curves and target height

bounds = function(theta)
{
  
  return(log_lik_lindley(theta, calltimes)-thres_height)
  
}

#Using uniroot() to find upper and lower bounds

theta_lower_bound = uniroot(bounds, interval= c(0.01, theta_predict))$root
theta_upper_bound = uniroot(bounds, interval = c(theta_predict, 1.0))$root

width_log= theta_upper_bound - theta_lower_bound

#defining the confidence interval

LR_CI = c(theta_lower_bound, theta_upper_bound)
#print(LR_CI)

cat("Log-Likelihood 95% Confidence Interval: [", LR_CI, "]")

Log-Likelihood 95% Confidence Interval: [ 0.1786338 0.265321 ]

theta_range = seq(0.1, 0.4, length.out = 100) #Data input range to cover the MLE

y_values = sapply(theta_range, bounds) #Calculate the 'height' for each theta using your bounds function #sapply tells R to run your function for every theta in the range


plot(theta_range, y_values, type = "l", lwd = 2, main = "LR Confidence Interval of Lindley Distrubution - Call times",
     xlab = "Theta",
     ylab = "")
abline(h=0, col = "blue")
theta_lower_bound = uniroot(bounds, interval= c(0.01, theta_predict))$root
theta_upper_bound = uniroot(bounds, interval = c(theta_predict, 1.0))$root

abline(v=c(theta_lower_bound, theta_upper_bound), col = "red")

Summarize results for question 1b \

The log likelihood distribution the telemarketer call time is estimated to be 0.2190994 and we are 95% confident the true call time is between [0.1786338, 0.265321]. The width between the confidence interval is 0.0866872. The interval includes the MLE Parameter 0.2190994, suggesting some accuracy for the confidence interval.

width_log < width_asym # FALSE log likelihood not a smaller range

[1] FALSE

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

The confidence intervals using both asymptotic and log likelihood are both comparable ranges that both hold the predicted MLE parameter 0.2190994 without , suggesting significance for both CI types. As the asymptotic confidence interval are [0.1758057, 0.2623931] with a width of [0.0865874] this is a slightly more narrow width than the log likelihood ratio CI at [0.1786338, 0.265321] with a width of [0.0866872].The slightly more narrow width of the asymptotic CI suggests the likelihood ratio CI cannot fully capture the true distribution shape of the Lindley distribution as well as the asymptotic CI can.

It should be noted that The log likelihood CI does not force shape as strictly as the asymptotic distribution that follow exact distributions such as standard normal, t or \(chi^2_1\) therefore log likelihood can better fit the natural distribution of the data set and are often more accurate than asymptotic normality, especially for a moderate sample size such as the one we have in this data set. In most cases the log likelihood CI ratio is more accurate than the asymptotic CI.

Considering the knowledge the log likelihood is usually the more accurate CI, the wider width in the log likelihood CI may suggest an balanced CI with more uncertainty rather than a lower quality in the CI itself. The overall recommendation is to use the log likelihood due to its resistance to being artificially narrow, unlike the asymptotic CI, and yield more accurate CI’s. Since the Lindley Distribution holds some skew in its distribution, using a CI type with the assumption normal or parametric shape, like the asymptotic CI, may lead to artificial CI widths due to its stronger reliance on approximation. The overall recommendation is to use the log-likelihood ratio CI for the call time data set due to its Lindley Distribution.