[R] parallel mle/optim and instability

From: Aaron J. Mackey <amackey_at_pcbi.upenn.edu>
Date: Thu 08 Jul 2004 - 22:58:57 EST

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 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Jul 08 23:06:44 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 02:35:15 EST