Hi

I am trying to recover the hessian of a problem optimised with
box-constraints. The problem is that in some cases, my estimates are very
close to the boundary, which will make optim(..., hessian=TRUE) or
optimHessian() fail, as they do not follow the box-constraints, and hence
estimate the function in the unfeasible parameter space.

As a simple example (my problem is more complex though, simple param
transformations do not apply ashere), imagine estimating mu and sigma
(restricted to be >0) of a simple normally distributed data, where however
the true sigma is very close to zero:

LLNorm <- function(para, dat=rn)  return(-sum(dnorm(dat, para[1], para[2],
log=TRUE)))
rn2 <- c(rep(10.3, 2000), 10.31)

>optim(c(10,1), fn=LLNorm, method="L-BFGS-B", lower=c(-Inf, 0.0000001),
dat=rn2,hessian=TRUE)
Error in optim(c(10, 1), fn = LLNorm, method = "L-BFGS-B", lower = c(-Inf,
: non-finite finite-difference value [2]

The only solution/workaround I found is to do a two steps procedure: use
optim() without hessian, then use optimHess, specifying the length of the
numerical variations (arg ndeps) as smaller as the distance between the
parameter and the bound, i.e.:
op<-optim(c(10,1), fn=LLNorm, method="L-BFGS-B", lower=c(-Inf, 0.0000001),
dat=rn2,hessian=FALSE)

H<-optimHess(op$par, LLNorm, dat=rn2,
control=list(ndeps=rep((op$par[2]-0)/1000,2)))

While this solution "works", it is surely not elegant, and I have following
questions:
1) is this solution to use smaller steps good at all?
2) I am not sure either I understand the reasons why the hessian in a
constrained optim problem is evaluated without the constraints, or at least
the problem transformed into a warning?

Furthermore, I realised that using the numDeriv package, the function
hessian() does not offer either constraints for the parameters, yet in the
special case of the example above, it does not fail (although would have
problems with parameter even closer to the bound), unlike the
optimHessian() function? See:
library(numDeriv)
hessian(LLNorm, op$par, dat=rn2)

This brings me to the final question: is there a technical reason not to
allow a "constrained" hessian, which seems to indicate the fact that the
two implementations in R do not do it? I found a very interesting answer by
Spencer Grave for  a similar question:
http://tolstoy.newcastle.edu.au/R/e4/help/08/06/15061.html
although this regards more the statistical implications than the numerical
issues.

Thanks a lot!

Matthieu

	[[alternative HTML version deleted]]

