#Define the negative log-likelihood function
nloglhood=function(p) 
  return( -(log(choose(100,10))+10*log(p)+90*log(1-p)) )
#Minimize the negative log-likelihood
binom.fit=optim(0.5,nloglhood,lower=0.0001,upper=0.9999,
                  hessian=T)
phat=binom.fit$par #The MLE

phat.var=1/binom.fit$hessian #Variance is inverse hessian
#Calculate approximate 95% Wald CI
phat+c(-1,1)*qnorm(0.975)*sqrt(phat.var)

library(Bhat) #Loading package Bhat
#Set up list for input into plkhci function
control.list=list(label="p",est=phat,low=0,upp=1)
#Calculate approximate 95% likelihood ratio CI
LRCI=plkhci(control.list,nloglhood,"p")

#Plot loglikelihood
p=seq(0.025,0.250,0.001)
loglhood=-nloglhood(p)
max.val=-binom.fit$value
plot(p,loglhood,,ylab="likelihood",xlim=c(0,max(p)),type="l",las=1)
segments(phat,-10,phat,max.val,lty=3)
#Add CI
crit.val=max.val-qchisq(0.95,1)/2
abline(h=crit.val,lty=2)
segments(LRCI[1],-10,LRCI[1],crit.val,lty=3)
segments(LRCI[2],-10,LRCI[2],crit.val,lty=3)
arrows(phat,crit.val,phat,max.val,length=0.10,code=3)
text(phat+0.01,max.val-1,"1.92",cex=0.8,srt=90)
arrows(LRCI[1],-8,LRCI[2],-8,length=0.10,code=3)
text(0.115,-7.8,"95% CI",cex=0.8)