#Negative log likelihood under H0: psi=psi0
nllhoodH0=function(lambda,y,psi0) {
  theta=c(psi0,lambda)
  return(nllhood(theta,y))
}
Faithful.fitH0=optim(c(80,5),nllhoodH0,y=waiting,psi0=c(1/3,55,5))
cat("Maximized log-lhood under H0:",-Faithful.fitH0$value,"\n")
cat("Unrestricted maximized log-lhood:",-Faithful.fit$value,"\n")
#Calculate likelihood ratio statistic
LRStat=2*(-Faithful.fit$value+Faithful.fitH0$value)
cat("LRT stat is",LRStat,"& p-value is",1-pchisq(LRStat,3),"\n")