lnLnorm <- function(theta,x,c) {
# compute lnL for normal distribution

# x is obs value or dl
# c is true if censored, F if observed
# theta is a vector of (mean, sd)

m <- theta[1]
s <- theta[2]
z <- (x-m)/s

lnl <- ifelse(c,log(pnorm(z)),-log(2*pi)/2 - log(s) -0.5*z^2)
sum(lnl)
}


zinc <- read.table('c:/philip/stat505/zinc.txt',as.is=T,header=T)

zinc.mle <- optim(c(9,4),lnLnorm,method='BFGS',hessian=T,control=list(fnscale=-1),
                  x=zinc$zinc,c=zinc$cens)
# arguments to optim are: starting values, 
#   function to be minimized
#   minimization method (BFGS is good general method, default is Nelder-Mead simplex
#   hessian: T if want to include the hessian
#   control parameters, fnscale=-1 multiplies function by -1
#     necessary because optim is a minimizer; we want a maximizer
#   x: vector of X values for lnLnorm
#   c: vector of C values for lnLnorm

# and here's the output:
#  $par is the mle's (mean, sd)
#  $value is the lnL at the mle
#  $counts is the number of times function was evaluated 
#  $convergence: 0 converged, 1 or other value = didn't
#  $hessian: numerical est. of 2nd deriv matrix               
zinc.mle
zinc.vc <- solve(-zinc.mle$hessian)
zinc.vc
          [,1]       [,2]
[1,]  0.9253705 -0.3712305
[2,] -0.3712305  0.8084565


# bad starting values can be bad
zinc.mle <- optim(c(1,1),lnLnorm,method='BFGS',hessian=T,control=list(fnscale=-1),
+                   x=zinc$zinc,c=zinc$cens)
>                   
> zinc.mle
$par
[1]  195.7277 2150.2950

$value
[1] -151.2024

$counts
function gradient 
     100      100 

$convergence
[1] 1

$message
NULL

$hessian
              [,1]         [,2]
[1,] -6.309619e-06 4.021672e-06
[2,]  4.021672e-06 3.232969e-06