[ADMB Users] The importance of parameterizing a nonlinear model well
dave fournier
davef at otter-rsch.com
Mon Aug 19 16:48:47 PDT 2013
I came across the following simple 3 parameter logistics problem
describing
how difficult it was to fit a small "weeds" data set.
https://groups.nceas.ucsb.edu/non-linear-modeling/projects/weeds/WRITEUP/weeds.pdf
there are 12 data points
t weeds
1 5.308
2 7.24
3 9.638
4 12.866
5 17.069
6 23.192
7 31.443
8 38.558
9 50.156
10 62.948
11 75.995
12 91.972
The model to be fit is
y = b_1/(1+b_2*exp(-b_3*t))
We are told several times how difficult this data set is to fit.
In fact is is a trivial data set to fit as I shall show. The
problem is that the parametrization used is a very poor one.
Consider the parameter b_1. This is the number of weeds at t=infinity.
It appears to be somewhat greater than 91, but how much? It is not
really obvious.
the parameter b_2 is even less intuitive.
we shall replace b_1 and b_2 with two parameters whose value is almost
obvious,
namely the number of weeds at time t_1=1, and t_n=12. Call thes
parameter y_1 and y_n.
clearly y_1 is close to 5.308 and y_n is close to 91.9
Given values for y_1, y_2, and b_3 it is a simple matter to solve for
b_2 and b_1.
The code in the tpl file is
b3=exp(logb3);
dvariable gamma1=exp(-b3*t(1));
dvariable gamman=exp(-b3*t(noObs));
b2=(yn-y1)/(y1*gamma1-yn*gamman); // solve for b2
b1=y1*(1.0+b2*gamma1); // solve for b1
For both models the sum of squares is rss 2.58728
Now lets make the data a bit more challenging. we change the 12'th data
point to
110.
For Ander's verion the sum of squares is rss 60.6799
while for my version it is rss 55.8948
Now change the 12'th data point to 125.
For Ander's version the sum of squares is rss 237.054
and for mine rss 133.582
So maybe we should think more about how to parameterize these models.
-------------- next part --------------
5.308 91.972 0 0
-------------- next part --------------
DATA_SECTION
init_int noObs
init_matrix obs(1,noObs,1,2)
PARAMETER_SECTION
init_number logb1;
init_number logb2;
init_number logb3;
init_number logSigma(-1);
sdreport_number sigma2;
sdreport_number b1;
sdreport_number b2;
sdreport_number b3;
sdreport_vector pred(1,noObs);
objective_function_value nll;
PRELIMINARY_CALCS_SECTION
PROCEDURE_SECTION
b1=exp(logb1);
b2=exp(logb2);
b3=exp(logb3);
sigma2=exp(2.0*logSigma);
for(int i=1; i<=noObs; ++i){
pred(i)=b1/(1.0+b2*exp(-b3*obs(i,1)));
// nll+=0.5*(log(2.0*M_PI*sigma2)+square((obs(i,2)-pred(i)))/sigma2);
}
nll+=noObs*0.5*log(norm2(column(obs,2)-pred));
REPORT_SECTION
report<<"rss " <<norm2(column(obs,2)-pred)<<endl;
GLOBALS_SECTION
#include <admodel.h>
void add_slave_suffix(const adstring& s){}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: weed.dat
Type: application/x-ns-proxy-autoconfig
Size: 214 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20130819/0c658e21/attachment.dat>
More information about the Users
mailing list