library(MASS)
data(Boston)
crim <- as.integer(Boston$crim > quantile(Boston$crim,0.8))
med.value <- Boston$medv
logitreg <- function(beta,x,y)
{   eta <- exp(beta[1] + beta[2]*x)
    p <- eta/(1+eta)
    return(-sum(ifelse(y,log(p),log(1-p)))) }
mle <- optim(c(0,0),logitreg,x=med.value,y=crim,hessian=TRUE)
se <- sqrt(diag(solve(mle$hessian)))
cbind(lower=mle$par - 1.96*se, upper=mle$par + 1.96*se)
a <- seq(-10,15,l=50)
b <- seq(-1,0.2,l=50)
L <- matrix(NA,50,50)
for(i in 1:50)
    for(j in 1:50)
        L[i,j] <- -logitreg(c(a[i],b[j]), med.value, crim)
par(mfrow=c(1,2))
contour(a,b,L)
persp(a,b,L)