[ADMB Users] The importance of parameterizing a nonlinear model well

[dave fournier]

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>