R-help,
I'm trying to estimate some parameters using the Maximum Likehood method.
The model describes fish growth using a sigmoidal-type of curve:
fn_w <- function(params) {
Winf <- params[1]
k <- params[2]
t0 <- params[3]
b <- params[4]
sigma <- params[5]
what <- Winf * (1-exp(- k *(tt - t0)))^b
logL <- -sum(dnorm(log(wobs),log(what),sqrt(sigma),TRUE))
return(logL)
}
tt <- 4:14
wobs <- c(1.545, 1.920, 2.321 ,2.591, 3.676, 4.425 ,5.028, 5.877, 6.990, 6.800
,6.900)
An then the optimization method:
OPT <-optim(c(8, .1, 0, 3, 1), fn_w, method="L-BFGS-B"
,lower=c(0.0, 0.001, 0.001,0.001, 0.01), upper = rep(Inf, 5), hessian=TRUE,
control=list(trace=1))
which gives:
$par Winf k t0 b
sigma
[1] 24.27206813 0.04679844 0.00100000 1.61760492 0.01000000
$value
[1] -11.69524
$counts
function gradient
143 143
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
$hessian
[,1] [,2] [,3] [,4] [,5]
[1,] 1.867150e+00 1.262763e+03 -7.857719 -5.153276e+01 -1.492850e-05
[2,] 1.262763e+03 8.608461e+05 -5512.469266 -3.562137e+04 9.693180e-05
[3,] -7.857719e+00 -5.512469e+03 41.670222 2.473167e+02 -5.356813e+01
[4,] -5.153276e+01 -3.562137e+04 247.316675 1.535086e+03 -1.464370e-03
[5,] -1.492850e-05 9.693180e-05 -53.568127 -1.464370e-03 1.730462e+04
after iteration number 80.
>From the biological point of view Winf =24(hipothesized asimptotical maximum
>weight)
makes no sense while the b parameter is no nearly close to b=3 leading to a
non-sigmoidal
curve.
However using the least-squares method provide with more sensible parameter
estimates
$par Winf k t0 b
[1] 10.3827256 0.0344187 3.1751376 2.9657368
$value
[1] 2.164403
$counts
function gradient
48 48
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
Is there anything wrong with my MLE function and parameters?
I have tried with distinct initial parameters also.
Can anyone give me a clue on this?
Thanks in advance.
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
