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

```