optim with contraints
On Saturday 21 June 2003 18:20, Spencer Graves wrote:
Adelchi: Permit me to add to Prof. Ripley's comments: If you want to know how to get around this kind of problem, I will tell you that I would modify the definition of "cp" to send the constraints to +/-Inf. Functions like "optim" tend to work better with unconstrained problems than constrained problems. In your example, I would typically replace the second element of the unknown in your example by its logarithm, then immediately exponentiate it as a first step of the function. Similarly I might replace the third by something like a logit transform z = log((1+x)/(1-x)), then immediately compute x = (1-exp(z))/(1+epx(z)) as a first step in the revised "cp". In many cases, these kinds of transformations have other advantages. Perhaps most importantly, a quadratic approximate tends to fit better over a wider range. This can reduce the number of iterations required for convergence. Moreover, if "cp" is a deviance = (-2)*log(likelihood), it can increase the accuracy of a normal approximation to the maximum likelihood estimates.
Spencer,
thanks for your comments. What you say is surely appropriate in most
circumstances. In fact, the original parametrisation of the problem was
unconstrained! and, yes, the function to be optimised is of type
deviance = (-2)*log(likelihood),
However, in this case, there are reasons to reparametrise in this form.
The essence is that near zero the unconstrained parameter has problems.
I can provide a detailed theoretical explanation of above statements,
but only if you really want it, and surely is not of interest to the whole
R-help list.
As Prof. Ripley indicated, R is a volunteer project. If you produce a reproducible example, it would help further if you step through the code for "optim" line by line and find where it bombs. If you can further suggest changes to the code that eliminate the problem, it could reduce substantially the time required for the bug to be fixed.
(On Monday) I shall see what I manage to understand: my R programming skill is
definetely not at "R developer" level. How do I step through the code for "optim"
line by line? its core statement is
res <- .Internal(optim(par, fn1, gr1, method, con, lower, upper))
and I do not know how to access the .Internal; probably these questions already
tell you a lot abut my skill..
I am aware of the fact that R is a volunteer project -- in fact a most appreciated
and fundamental project of the statistical community. In writing my earlier
message, I did intended to imply that some member of R-core had to rush looking
at my problem.
With best wishes,
Adelchi
Adelchi Azzalini <azzalini at stat.unipd.it> Dipart.Scienze Statistiche, Universit? di Padova, Italia http://azzalini.stat.unipd.it/