[R] optimization problem

[Hector Villalobos]

Dear R users,

I'im trying to find the parameters of a dynamic biomass model using maximum
likelihood estimation. I used two approaches, one by hand, with optim()
function and the other using mle2() function from package bbmle. My problem
is that the results change a lot depending on the initial values... I can't
see what I am doing wrong...

Thank you for any help!


# Data
x <- 1995:2010
B <- c(3500, 3200, 3000, 2800, 2600, 3000, 3200, 3800, 4200, 4300, 4400,
4400, 4500, 4600, 5000, 4300)
Ct <- c(912, 767, 642, 482, 353, 331, 332, 309, 366, 402, 392, 478, 408,
434, 407, 637)
a <- c(0.539, 0.603, -0.948, 0.166, 1.895, 0.786, 0.901, 0.844, 0.337,
0.429, 0.304, 0.230, 1.001, 0.750, 0.507, 1.502)
Ag <- 0.55

# Function with quantity to minimize
modl <- function(par) {
  ro <- par[1]
  ko <- par[2]
  n <- length(B)
  Be <- rep(NA, n)
  Be[1] <- ko * Ag
  for ( k in 2:n)
    Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]
  err <- (log(B) - log(Be))^2
  ee <- sqrt( sum(err)/(n-2) )
  LL <- (1/(sqrt(2*pi)*ee)) * exp( -(err/(2*ee^2) ) )
  -crossprod(LL)
}

# Using function optim()
par.optim <- optim(par = list(ro=0.4, ko=8000), modl, method = "BFGS")
ro <- par.optim$par[1]
ko <- par.optim$par[2]

# estimated values of "B"
n <- length(B)
Be <- rep(NA, n)
Be[1] <- ko * Ag
for ( k in 2:n)
  Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]

# Plot, estimation of "B" seems reasonable....
plot(x, B, ylim=c(1000, 7000))
lines(x, Be, col="blue", lwd=2)


# ... but it is very sensible to initial values...
par.optim2 <- optim(par = list(ro=0.4, ko=10000), modl, method = "BFGS")
ro2 <- par.optim2$par[1]
ko2 <- par.optim2$par[2]

Be2 <- rep(NA, n)
Be2[1] <- ko2 * Ag
for ( k in 2:n)
  Be2[k] <- Be2[k-1] + ro2 * a[k-1] * Be2[k-1] * (1 - Be2[k-1]/ko2) -
Ct[k-1]

lines(x, Be2, col="blue", lwd=2, lty=3)



# Uing mle2 function
library(bbmle)
LL <- function(ro, ko, mu, sigma) {
  n <- length(B)
  Be <- rep(NA, n)
  Be[1] <- ko * Ag
  for ( k in 2:n)
    Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]
  err <- log(B) - log(Be)
  R <- (dnorm(err, mu, sigma, log=TRUE))
  -sum(R)
}

Bc.mle <- mle2(LL, start = list(ro=0.4, ko=8000, mu=0, sigma=1))
summary(Bc.mle)

ro3 <- coef(Bc.mle)[1]
ko3 <- coef(Bc.mle)[2]

Be3 <- rep(NA, n)
Be3[1] <- ko3 * Ag
for ( k in 2:n)
  Be3[k] <- Be3[k-1] + ro3 * a[k-1] * Be3[k-1] * (1 - Be3[k-1]/ko3) -
Ct[k-1]

lines(x, Be3, col="red", lwd=2)


-- 

Héctor Villalobos <Hector.Villalobos.O at gmail.com>

Depto. de Pesquerías y Biología Marina

Centro Interdisciplinario de Ciencias Marinas-Instituto Politécnico Nacional

CICIMAR-I.P.N.

A.P. 592. Colonia Centro

La Paz, Baja California Sur, MÉXICO. 23000

Tels.: (+52 612) 122 53 44; 123 46 58; 123 47 34 ext.: 81602

Fax: (+52 612) 122 53 22

	[[alternative HTML version deleted]]