[R] parallel mle/optim and instability
Spencer Graves
spencer.graves at pdf.com
Thu Jul 8 18:17:00 CEST 2004
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
More information about the R-help
mailing list