[R] parallel mle/optim and instability

Aaron J. Mackey amackey at pcbi.upenn.edu
Thu Jul 8 14:58:57 CEST 2004


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",
          );




More information about the R-help mailing list