Re: [R] optim() not finding optimal values

From: dave fournier <otter_at_otter-rsch.com>
Date: Mon, 28 Jun 2010 04:09:49 -0400

If you are going to make this program available for general   use you want to take every precaution to make it bulletproof.

This is a fairly informative data set. The model will undoubtedly   be used on far less informative data. While the model looks   pretty simple it is very challenging from a numerical point of view.   I took a moment to code it up in AD Model Builder. The true minimum is   1619.480495 So I think Ravi has finally arrived pretty close to the answer.
One way of judging the difficulty of a model is to look at the   eigenvalues of the Hessian at the minimum. They are

3.122884668e-09 1.410866202e-08 1866282.520 1.330233652e+13

so the condition number is around 1.e+21. One begins to see why these   models are challenging. The model as formulated represents the state   of the art in fisheries models circa 1985.   A lot of progress has been made since that time.   Using B_t for the biomass and C_t for the catch the equation   in the code is

B_{t+1} = B_t + r *B_t*(1-B_t/K) -C_t (1)   First notice that
for (1) to make sense the following conditions must be satisfied

B_t > 0 for all t
r > 0
K>0

Strictly speaking it is not necessary that B_t<=K but if B_t>K and r   is large then B_{t+1} could be <0. So formulation (1) gives   Murphys law a good chance. How to fix it. Notice that (1) is really   a rough approximation to the solution of a differential equation

B'(t) = r *B(t)*(1-B(t)/K) -C (2)

where in (2) C is a constant catch rate. To fix (1) we use   a semi-implicit differencing scheme. Because it is useful to   allow smaller step sizes than one we denote them by d.

B_{t+d} = B_t + d* r *B_t*(1-B_{t+d}/K) -d*C_t*B_{t+d}/B_t (1)

The idea is that the quantity 1-x with x>0 will be replaced by   1/(1+x). Expanding 2 and solving for B_{t+d} yields

B_{t+d} = (1+d*r) B_t / (1+d*r*B_t/K +d*C_t/B_t) (3)

So long as r>0, K>0 C_t>0 then starting from an initial value    B_0 > 0 ensures that B_t> 0 for all t>0. We can let    d=1/nsteps where nsteps is the number of steps in the    approximate integration for each year    which can be increased until the solution is judged to be close    enough to the exact solution from (2)

Notice that in (3) as C_t --> infinity B_{t+d} --> 0    So that you can never catch more fish than you have.

I coded up this version of the model in AD Model Builder and    fit it to the data. It is now much more resistant to bad    starting values for the parameters etc.

If anyone wants the tpl file for the model in ADMB they can    contact me off list.

R-help_at_r-project.org 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. Received on Mon 28 Jun 2010 - 16:02:46 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 28 Jun 2010 - 16:30:43 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.