#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
plkhci(control.list,nloglhood,"p")