[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