[R] optim with contraints

[Spencer Graves]

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.

	  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.

hth, spencer graves

Prof Brian Ripley wrote:
> Adelchi,
> 
> R is a volunteer project, and we would need a volunteer to look at this.  
> `The developers of optim' (who are R-core) used published code for this
> method. Since you are getting different answers on different platforms
> this might be a bug in your compiler or run-time rather than R.
> 
> Please submit a reproducible bug report to R-bugs.  As you will see some
> bugs don't get fixed very fast (including one I submitted from Padua in
> 2001).
> 
> Brian
> 
> On Sat, 21 Jun 2003, Adelchi Azzalini wrote:
> 
> 
>>There seems to exist peculiar cases where optim does not take care
>>of constraints on the parameters to be optimized over.  The call to
>>optim is of the form
>>
>>  opt <- optim(cp, fn=sn.dev, gr=sn.dev.gh, method="L-BFGS-B",
>>           lower=c(-Inf, 1e-10, -0.99527), 
>>           upper=c( Inf, Inf,    0.99527), 
>>           control=control, X=X, y=y, hessian=FALSE)
>>
>>The code has worked fine many times, but I have come across cases 
>>(for suitable data X and y) where the constraint on the last component 
>>is ignored;  that means that a call is made to sn.dev with
>>       cp = -1.3546  0.4645  3.1741
>>so the third component exceeds 0.99527, and the program stops 
>>because of a check in the  function to be obtimised.
>>
>>The call just before was to the gradient function
>>sn.dev.gh: gradient =    219013   -312643 441647332
>>which has rather large values.
>>
>>To make the problem more interesting, it shows up on a Linux
>>(Debian) installation, but it works fine on MS-windows.
>>In both cases, R is 1.7.0.  
>>
>>Perhaps this sort of question should not be directed to the R-help
>>list, rather to the developer(s) of optim. Please instruct me on this 
>>point. Also, I appreciate that the above is not a reproducible example;
>>this would be longish in text, and I though to ask first to whom 
>>it is appropriate that I direct my question.
>>
>>
>>Best wishes,
>>
>>Adelchi Azzalini
>>
>>
>>
> 
>