Your model is producing -Inf entries in the vector Be (in function modl and LL)
at some stage during the optimization process.
You should first do something about that before anything else.
Berend
> On 17 Oct 2015, at 03:01, Bert Gunter <bgunter.4...@gmail.com> wrote:
>
> I made no attempt to examine your details for problems, but in general,
>
> My problem
>> is that the results change a lot depending on the initial values... I can't
>> see what I am doing wrong...
>>
>> This is a symptom of an overparameterized model: The parameter estimates
>> are unstable even though the predictions may not change much. In other
>> words, your model may be too complex for your data.
>
>
> Whether that is true here, you or others will have to determine. Try
> simplifying your model as a start.
>
> -- Bert
>
>
>
>>
>>
>> # 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.villalobo...@gmail.com <javascript:;>>
>>
>> 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]]
>>
>> ______________________________________________
>> R-help@r-project.org <javascript:;> mailing list -- To UNSUBSCRIBE and
>> more, see
>> 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.
>
>
>
> --
> Bert Gunter
>
> "Data is not information. Information is not knowledge. And knowledge is
> certainly not wisdom."
> -- Clifford Stoll
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
