[R] problem with convergence in mle2/optim function
Adam Zeilinger
zeil0006 at umn.edu
Fri Oct 5 07:12:07 CEST 2012
Hello R Help,
I am trying solve an MLE convergence problem: I would like to estimate
four parameters, p1, p2, mu1, mu2, which relate to the probabilities,
P1, P2, P3, of a multinomial (trinomial) distribution. I am using the
mle2() function and feeding it a time series dataset composed of four
columns: time point, number of successes in category 1, number of
successes in category 2, and number of success in category 3. The
column headers are: t, n1, n2, and n3.
The mle2() function converges occasionally, and I need to improve the
rate of convergence when used in a stochastic simulation, with multiple
stochastically generated datasets. When mle2() does not converge, it
returns an error: "Error in optim(par = c(2, 2, 0.001, 0.001), fn =
function (p) : L-BFGS-B needs finite values of 'fn'." I am using the
L-BFGS-B optimization method with a lower box constraint of zero for all
four parameters. While I do not know any theoretical upper limit(s) to
the parameter values, I have not seen any parameter estimates above 2
when using empirical data. It seems that when I start with certain
'true' parameter values, the rate of convergence is quite high, whereas
other "true" parameter values are very difficult to estimate. For
example, the true parameter values p1 = 2, p2 = 2, mu1 = 0.001, mu2 =
0.001 causes convergence problems, but the parameter values p1 = 0.3, p2
= 0.3, mu1 = 0.08, mu2 = 0.08 lead to high convergence rate. I've
chosen these two sets of values because they represent the upper and
lower estimates of parameter values derived from graphical methods.
First, do you have any suggestions on how to improve the rate of
convergence and avoid the "finite values of 'fn'" error? Perhaps it has
to do with the true parameter values being so close to the boundary? If
so, any suggestions on how to estimate parameter values that are near zero?
Here is reproducible and relevant code from my stochastic simulation:
########################################################################
library(bbmle)
library(combinat)
# define multinomial distribution
dmnom2 <- function(x,prob,log=FALSE) {
r <- lgamma(sum(x) + 1) + sum(x * log(prob) - lgamma(x + 1))
if (log) r else exp(r)
}
# vector of time points
tv <- 1:20
# Negative log likelihood function
NLL.func <- function(p1, p2, mu1, mu2, y){
t <- y$tv
n1 <- y$n1
n2 <- y$n2
n3 <- y$n3
P1 <- (p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu2*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu2*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))))
P2 <- (p2*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*(-mu1^2 + 2*mu2*p1 +
mu1*(mu2 - p1 + p2)) - mu1*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu1*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu1*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))))
P3 <- 1 - P1 - P2
p.all <- c(P1, P2, P3)
-sum(dmnom2(c(n1, n2, n3), prob = p.all, log = TRUE))
}
## Generate simulated data
# Model equations as expressions,
P1 <- expression((p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu2*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu2*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))
P2 <- expression((p2*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*(-mu1^2 + 2*mu2*p1 +
mu1*(mu2 - p1 + p2)) - mu1*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu1*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu1*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))
# True parameter values
p1t = 2; p2t = 2; mu1t = 0.001; mu2t = 0.001
# Function to calculate probabilities from 'true' parameter values
psim <- function(x){
params <- list(p1 = p1t, p2 = p2t, mu1 = mu1t, mu2 = mu2t, t = x)
eval.P1 <- eval(P1, params)
eval.P2 <- eval(P2, params)
P3 <- 1 - eval.P1 - eval.P2
c(x, matrix(c(eval.P1, eval.P2, P3), ncol = 3))
}
pdat <- sapply(tv, psim, simplify = TRUE)
Pdat <- as.data.frame(t(pdat))
names(Pdat) <- c("time", "P1", "P2", "P3")
# Generate simulated data set from probabilities
n = rep(20, length(tv))
p = as.matrix(Pdat[,2:4])
y <- as.data.frame(rmultinomial(n,p))
yt <- cbind(tv, y)
names(yt) <- c("tv", "n1", "n2", "n3")
# mle2 call
mle.fit <- mle2(NLL.func, data = list(y = yt),
start = list(p1 = p1t, p2 = p2t, mu1 = mu1t, mu2 = mu2t),
control = list(maxit = 5000, factr = 1e-10, lmm = 17),
method = "L-BFGS-B", skip.hessian = TRUE,
lower = list(p1 = 0, p2 = 0, mu1 = 0, mu2 = 0))
###########################################################################
I interpret the error as having to do with the finite difference
approximation failing. If so, perhaps a gradient function would help?
If you agree, I've described my unsuccessful attempt at writing a
gradient function below. If a gradient function is unnecessary, ignore
the remainder of this message.
My gradient function: I derived the gradient function by taking the
derivative of my NLL equation with respect to each parameter. My NLL
equation is the probability mass function of the trinomial distribution.
Thus the gradient equation for, say, parameter p1 would be:
gr.p1 <- deriv(log(P1^n1), p1) + deriv(log(P2^n2), p1) +
deriv(log(P3^n3), p1)
This produces a very large equation, which I won't reproduce here.
Let's say that the four gradient equations for the four parameters are
defined as gr.p1, gr.p2, gr.mu1, gr.mu2, and all are derived as
described above for gr.p1. These gradient equations are functions of
p1, p2, mu1, mu2, t, n1, n2, and n3. My current gradient function is:
grr <- function(p1, p2, mu1, mu2, y){
t <- y[,1]
n1 <- y[,2]
n2 <- y[,3]
n3 <- y[,4]
gr.p1 <- .......
gr.p2 <- .......
gr.mu1 <- .......
gr.mu2 <- .......
c(gr.p1, gr.p2, gr.mu1, gr.mu2)
}
The problem is that I need to supply values for t, n1, n2, and n3 to the
gradient function, which are from the dataset yt, above. When I supply
the dataset yt, the function produces a vector of length 4*nrow(yt) =
80. When I include it in my mle2() function, I get an error that mle2
(optim) requires a vector of length 4. How do I write my gradient
function to work in mle2()?
Any help would be much appreciated.
Adam Zeilinger
--
Adam Zeilinger
Post Doctoral Scholar
Department of Entomology
University of California Riverside
www.linkedin.com/in/adamzeilinger
More information about the R-help
mailing list