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