Thanks to all who responded, I've found a very useful code here:
http://courses.washington.edu/fish507/notes.html In particular the Lecture 3... Héctor 2015-10-17 7:05 GMT+00:00 Berend Hasselman <b...@xs4all.nl>: > > 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. > > -- Dr. Héctor Villalobos <hector.villalobo...@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]] ______________________________________________ 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.
