est.bes <-
function(Y,X,intercept=TRUE,start.vals=rep(0,ncol(X)+as.numeric(intercept))){

x<-start.vals
opt<-optim(x,fn=lr.loglike,gr=NULL,Y=Y,X=X,intercept=intercept,method="BFGS",control=list(fnscale=-1),hessian=T)
beta.hat<-opt$par
beta.cov<-solve(opt$hessian)/length(Y)
if(intercept){
Xm<-cbind(rep(1,length(Y)),X)
}else{
Xm<-X
}
fitted.lp<-Xm %*% beta.hat
fitted.Y<-expit(fitted.lp)
conv<-TRUE
if(opt$convergence!=0){
conv<-FALSE
print('Optim did not converge')
}

return(list(ests=beta.hat,varcov=beta.cov,fitted.lp=fitted.lp,fitted.Y=fitted.Y,conv=conv))
}