#This is R code that uses observed exponential data to find a 2-sided,
#90% CI by inverting the LRT. You should run this code line-by-line and
#interactively to see how powerful a tool R is for understanding concepts.
#The code works from example 9.2.3 in the text to obtain the CI.
#However, unlike the example in the text, it implements empirical
#approaches to finding both the critical value k of the test and the CI.
#Here is the real data we will use, assumed to come from an exponential distribution.
x=c(1.47,1.84,1.49,3.62,.86,1.02,2.65,1.55,13.84,.52)
######################################################################
#STEP 1: Get the critical value k for the acceptance region of a
#level-alpha LRT. This step does not depend on the observed data x.
alpha=.10 #The level alpha of the test to be inverted.
#Under H0 that lambda=lambda0, we know that Y=sum(X1,..,X10)/lambda0
#is a gamma rv with shape parameter n and scale parameter 1.
#To find the value of k, we'll take an empirical approach as follows.
#Use R to generate a huge number of realizations (one million) of the rv Y.
largeN=1000000
yscale=1
n=length(x)
help(rgamma)
Y=rgamma(largeN, shape=n, rate=1/yscale)
#Let W=(Y^n)*exp(-Y) to mimic the left-hand side of the text's
#equation 9.2.2 for the acceptance region
W=Y^n*exp(-Y)
#Plot the null distribution of W and look at it.
hist(W,breaks="FD", main= "Null density for determining k")
#Want to find a k such that the probability mass of W to the right
#of k is 1-alpha under the null distribution. Thus k is the alpha-th
#quantile of the null distribution of W.
help(quantile)
k=quantile(W,alpha)
help(abline)
#Use a vertical line to superpose k on the plot of the distribution of W
abline(v=k,lty=2)
#Should check that your value of k agrees with the value obtained by others
k
## 10%
## 114671.7
##############################################################
#STEP 2: Now, using the critical value k obtained in the first
#step, and the observed data x, find all lambda0 values such that
#the observed data accepts H0:lambda=lambda0
#Make candidate lambda0's run from near 0 to near 20 on a fine grid
lambda0=(1:20000)/1001 #there are 20000 candidate lambda0's
#Consider the event that the data X are in the acceptance region of
#the level-alpha LRT. Get the left-hand-side of the expression that
#defines this event in equation 9.2.2.
lefthandside=((sum(x)/lambda0)^n ) * exp(- sum(x)/lambda0)
#Check if lefthandside is greater than k for any of the lambda0 values tried
accept=(lefthandside>k)
#From the set of all candidate lambda0's, pluck out the ones that
#the observed data x accepted. The range of these accepted lambda0's
#is the 1-alpha * 100% CI
CI=range(lambda0[accept])
#Check what is the value of the empirically-determined CI
CI
## [1] 1.782218 5.126873