From: dave fournier <otter_at_otter-rsch.com>

Date: Mon, 28 Jun 2010 04:09:49 -0400

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

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

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.
*