[R] optim with contraints

[Adelchi Azzalini]

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/