[R] parallel mle/optim and instability

[Spencer Graves]

      Does "optim" give you "computationally singular" errors?  I've had 
similar problems with "nls", but "optim" have given me answers even in 
such cases. 

      Do your numbers have substantially different orders of magnitude?  
Is it feasible to rescale everything to mean 0, standard deviation of 1, 
and then write your model in terms of these rescaled variables?  If yes, 
this can help. 

      I suggest you reparameterize the problem to move the constraints 
to +/-Inf.  If you can't do that, add penalty function corresponding to 
the constraints, while making sure your objective function will never 
give NA nor +/-Inf.  Find the conditions that generate that and 
eliminate them. 

      Without box constraints, "optim" has 4 methods.  I'd try the first 
3;  the fourth is simulated annealing, which may be wonderful for 
objective functions with many local minima but otherwise takes forever 
and produces answers that are inferior to the other three methods, in 
the few cases I've tried. 

      Also, with a problem with many parameters, I've written a function 
that allows me to pass as a "..." argument the subset of parameters I 
want to estimate.  Then I can try different subsets until I isolate the 
problem. 

      hope this helps.  spencer graves

Aaron J. Mackey wrote:

>
> I have a MLE task that for a small number of parameters finishes in a 
> reasonable amount of time, but for my "real" case (with 17 parameters 
> to be estimated) either takes far too long (over a day), or fails with 
> "computationally singular" errors.  So a) are there any parallel 
> implementations of optim() (in R or otherwise) and b) how can I make 
> my function more robust? (I've already resorted to using bounded 
> parameters and log transformations, which has helped to some degree)
>
> Thanks, -Aaron
>
> library(stats4);
> d <- data.frame( ix=c(0,1,2,3,4,5,6,7),
>                  ct=c(1000,9609,18403,2617,8237,3619,5520,18908),
>                  A=c(0,1,0,1,0,1,0,1),
>                  B=c(0,0,1,1,0,0,1,1),
>                  C=c(0,0,0,0,1,1,1,1)
>                );
> ct <- round(logb(length(d$ix), 2))
> ll <- function( th=0.5,
>                 a1=log(0.5), a2=log(0.5), a3=log(0.5),
>                 b1=log(0.5), b2=log(0.5), b3=log(0.5)
>               ) {
>   a <- exp(sapply(1:ct, function (x) { get(paste("a", x, sep="")) }));
>   b <- exp(sapply(1:ct, function (x) { get(paste("b", x, sep="")) }));
>   s <- -sum( d$ct * log( sapply( d$ix,
>                                  function (ix, th, a, b) {
>                                     x <- d[ix+1,3:(ct+2)]
>                                     (th     * prod((b ^ (1-x)) * 
> ((1-b) ^ x    ))) +
>                                     ((1-th) * prod((a ^ x    ) * 
> ((1-a) ^ (1-x))))
>                                   },
>                                   th, a, b
>                                )
>                        )
>         );
>   if (!is.finite(s)) stop(cat(th, a, b, s, sep="\n"));
>   s;
> }
>
> ml <- mle(ll,
>           lower=c(    1e-10, rep(log(    1e-10), 2*ct)),
>           upper=c(1 - 1e-10, rep(log(1 - 1e-10), 2*ct)),
>           method="L-BFGS-B",
>          );
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html