[R] optimization problem

[Hector Villalobos]

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 <bhh at 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.4567 at 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.Villalobos.O at 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 at 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 at 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.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]]