#Draws a likelihood ratio CI plot
#Required arguments are
# neg log likelihood function
# MLE (i.e., theta hat)
# upper and lower bounds of LRCI
# sequence of theta values over which to evaluate the log-likelihood
LRCIplot=function(nllhood,MLE,LRCI,thseq,xlab=NULL,ylab=NULL,...) {
  if(is.null(xlab)) xlab=expression(theta)
  if(is.null(ylab)) ylab="Log-likelihood"
  lhoodseq=numeric(length(thseq))
  for(i in 1:length(thseq)) lhoodseq[i]=-nllhood(thseq[i],...)
  plot(thseq,lhoodseq,type="l",las=1,xlab=xlab,ylab=ylab)
  critval=-nllhood(MLE)-qchisq(0.95,1)/2
  abline(h=critval,lty=2)
  segments(LRCI[1],min(lhoodseq)-10,LRCI[1],critval,lty=3)
  segments(LRCI[2],min(lhoodseq)-10,LRCI[2],critval,lty=3)
  arrows(MLE,critval,MLE,-nllhood(MLE),length=0.10,code=3)
  text(MLE-0.02*(max(thseq)-min(thseq)),-nllhood(MLE)-1,"1.92",cex=0.8,srt=90)
  arrow.y.posn=min(lhoodseq)+0.2*(critval-min(lhoodseq))  
  arrows(LRCI[1],arrow.y.posn,LRCI[2],arrow.y.posn,length=0.10,code=3)
  text.y.posn=arrow.y.posn+0.03*(max(lhoodseq)-min(lhoodseq))
  text(mean(LRCI),text.y.posn,"95% CI",cex=0.8) }