###################################################
#
# Bootstrap for the Censored Poisson Example
# 
# Here, verifying the distribution of 2[log(LR)]
# against asymptotic theory.
#
# Professor M. Zhu
# Current version  = Nov 2022
#
###################################################

rm(list=ls(all=TRUE))

myloglik <- function(theta, v, y)
{
# computes the log-likelihood function for the
# censored Poisson model, given data (v,y) ...
# written by M. Zhu for teaching illustration
  u <- exp(-v*theta)
  return(sum(y*log(1-u)-(1-y)*v*theta))
}

# Here are some data:
v<-c(8,8,8,8,8,8,8,8,8,8,
     4,4,4,4,4,4,4,4,4,4,
     2,2,2,2,2,2,2,2,2,2,
     1,1,1,1,1,1,1,1,1,1)
y<-c(1,1,1,1,1,1,1,1,1,1,
     1,1,1,1,1,1,1,1,0,0,
     1,1,1,1,1,1,1,0,0,0,
     1,1,1,0,0,0,0,0,0,0)

# data estimate (here, using a built-in optimizer, 
# rather than writing our own Newton's procedure, but
# can certainly use our own Newton's procedure as well)
thetahat<-optimize(myloglik,interval=c(0,1),maximum=TRUE,
                                             v=v,y=y)$maximum
#######################
# do the bootstrap here
#######################
B=1000
logLRstar<-thetastar<-numeric(B)
for (b in 1:B){
 ind=sample(1:40,size=40,replace=TRUE)
 thetastar[b]<-optimize(myloglik,interval=c(0,1),maximum=TRUE,
                                             v=v[ind],y=y[ind])$maximum
 logLRstar[b]=2*(myloglik(theta=thetastar[b],v=v[ind],y=y[ind])
                -myloglik(theta=thetahat,    v=v[ind],y=y[ind]))
}

##############################
# code below for making a plot
##############################

# set graphic parameter limits
hi=ceiling(max(logLRstar))

# histogram of bootstrap samples
hist(logLRstar,xlim=c(0,hi),freq=FALSE,xlab='',ylab='',main='',
  border=FALSE,nclass=15)->junk

# add theoretical distribution based on asymptotic arguments, chi-squared w/ df=1
xx=seq(0,hi,len=100)
lines(xx,dchisq(xx,df=1),col='red',lwd=2)

# add quantiles and text label
abline(v=quantile(logLRstar,prob=0.95),col='black',lty=2)
abline(v=qchisq(0.95,df=1),col='red',lty=2)
text(4,max(junk$density)-0.025,'95% quantiles',pos=4,col='blue')

# add title to graphics
title(main='Bootstrap vs Asymptotic Distributions of 2[log(LR)]')