# Re: [R] parallel mle/optim and instability

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Fri 09 Jul 2004 - 02:17:00 EST

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@stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help