n<-1000
set.seed(2342)

X<-rnorm(n,0,2)
be.vec<-c(-0.5,1.0)

eta<-cbind(rep(1,n),X) %*% be.vec

pvec<-1/(1+exp(-eta))
Y<-rbinom(n,1,pvec)

table(Y)

log.like<-function(bev,Xv,Yv){

	ev<-cbind(rep(1,length(Yv)),Xv) %*% bev
	pv<-1/(1+exp(-ev))
	tvec<-sum(dbinom(Yv,1,pv,log=T))
	return(tvec)

}



log.like(be.vec,X,Y)
be0<-be1<-seq(-1,3,by=0.1)
lmat<-matrix(0,nrow=length(be0),ncol=length(be1))
for(i in 1:length(be0)){
	for(j in 1:length(be1)){
	
		lmat[i,j]<-log.like(c(be0[i],be1[j]),X,Y)

	}
	
}

image(be0,be1,lmat,col=heat.colors(50))
contour(be0,be1,lmat,nlevels=50,add=T)
abline(v=be.vec[1],h=be.vec[2],col="red")

image(be0,be1,lmat,col=heat.colors(50),xlab=expression(beta[0]),ylab=expression(beta[1]))
contour(be0,be1,lmat,nlevels=50,add=T)
abline(v=be.vec[1],h=be.vec[2],col="red")

######################################################################
optim(c(0,0),fn=log.like,Yv=Y,Xv=X,control=list(fnscale=-1))

logreg.opt<-optim(c(0,0),fn=log.like,Yv=Y,Xv=X,
                  method="BFGS",control=list(fnscale=-1),hessian=T)


Asymp.var<-solve(-logreg.opt$hessian)
                         


library(MASS)
Approx.Beta<-mvrnorm(2000,mu=logreg.opt$par,Sigma=Asymp.var)


be0.new<-seq(-1,0,by=0.01)
be1.new<-seq(0.5,1.5,by=0.01)

lmat2<-matrix(0,nrow=length(be1.new),ncol=length(be0.new))
for(i in 1:n){
	lpmat<-outer(be0.new,be1.new*X[i],"+")
	lmat2<-lmat2+Y[i]*lpmat-log(1+exp(lpmat))
}

image(be0.new,be1.new,lmat2,col=heat.colors(50),xlab=expression(beta[0]),ylab=expression(beta[1]))
contour(be0.new,be1.new,lmat2,add=T,nlevels=50)
points(Approx.Beta,pch=19,cex=0.5)