simData = function(X,beta) {
   "simuate y~Bernoulli(X beta)"
   eta = X %*% matrix(beta,ncol=1)
   pvec = 1.0/(1.0 + exp(-eta))
   y = rbinom(nrow(X),1,pvec)
   return(y)
}

phat = function(X,beta) {
   "p(X beta)"
   eta = X %*% matrix(beta,ncol=1)
   pvec = 1.0/(1.0 + exp(-eta))
   return(pvec)
}

mLL = function(X,y,beta) {
   eta = X %*% matrix(beta,ncol=1)
   py1 = 1.0/(1.0 + exp(-eta))
   return(-sum(ifelse(y==1,log(py1),log(1-py1)))/n)
}

lgGrad = function(X,y,beta) {
   eta = X %*% matrix(beta,ncol=1)
   pvec = 1.0/(1.0 + exp(-eta))
   gvec = t(X) %*% matrix(y-pvec,ncol=1)
   return (-gvec/n)
}

lgH = function(X,beta) {
   eta = X %*% matrix(beta,ncol=1)
   pvec = as.double(1.0/(1.0 + exp(-eta)))
   D = diag(pvec*(1-pvec))
   Hm = (t(X) %*% D %*% X)/n
   return(Hm)
}