In the course of revising a paper I have had occasion to attempt to
maximize a rather
complicated log likelihood using the function nlm(). This is at the
demand of a referee
who claims that this will work better than my proposed use of a home-
grown implementation
of the Levenberg-Marquardt algorithm.
I have run into serious hiccups in attempting to apply nlm().
If I provide gradient and hessian attributes to the returned function
value
(the ***negative*** of the log likelihood, since nlm() minimizes
things) then
nlm() wanders around for a very long time and reaches a highly sub-
optimal
value of the negative log likelihood.
It also gives multiple warnings to the effect ``NA/Inf replaced
by maximum positive value in nlm(etc.)''. Putting a browser() inside
my objective
function appears to reveal that the hessian becomes singular in many
calls to
the function. I presume (presumptuously) that this is the cause of
the NA/Inf-s.
Also if check.analyticals is TRUE I get an error to the effect
``probable coding
error in analytic hessian''.
Since the coding of the hessian appears to work ``perfectly'' with my
home-grown
Levenberg-Marquardt procedure (I get answers completely consistent
with those
obtained from other approaches including use of the EM algorithm and
use of
numerical optimization invoking optim()) I am at least 99% confident
that my coding
of the hessian is NOT incorrect.
If I do not provide gradient and hessian attributes but simply let nlm
() do its thing
completely numerically, then the algorithm converges to (within
numerical
tolerance) the same answer obtained using my other three methods (i.e.
Levenberg-Marquardt, EM, and optim()). However I still get a warning
(one only)
in respect of NA/Inf being replaced by maximum positive value.
The hessian returned by nlm() in this case agrees closely with that
calculated
by my ``analytical'' procedure, which gives me further confidence
that my calculation
of the hessian is not incorrect despite the claim induced by
``check.analyticals=TRUE''.
My gut feeling is that the problem is just a bit too complicated for
nlm() to handle, and
I would simply leave it and go with the other approaches which appear
to work consistently
and well. However Referees Must Be Satisfied. Has anyone any
suggestions as to
how I can get nlm() to work with the analytical gradient and hessian?
The code is much too complicated to include in this posting, but if
anyone is interested
in experimenting I can make available a package that I have prepared
which does
all of the calculations in a convenient manner and demonstrates the
problem.
Thanks for any help that anyone can give.
cheers,
Rolf Turner
######################################################################
Attention:\ This e-mail message is privileged and confidenti...{{dropped}}