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))
}