$ontext

  Nonlinear Least Squares Regression example

  Erwin Kalvelagen, nov 2007

  Reference:
       http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml

--------------------------------------------------------------------------
Procedure:     Nonlinear Least Squares Regression

Description:   This model and data are an example of fitting
               sigmoidal growth curves taken from Ratkowsky (1983).
               The response variable is pasture yield, and the
               predictor variable is growing time.

Reference:     Ratkowsky, D.A. (1983).
               Nonlinear Regression Modeling.
               New York, NY:  Marcel Dekker, pp. 61 and 88.

Data:          1 Response  (y = pasture yield)
               1 Predictor (x = growing time)
               9 Observations
               Higher Level of Difficulty
               Observed Data

Model:         Exponential Class
               3 Parameters (b1 to b3)

               y = b1 / (1+exp[b2-b3*x])  +  e

          Starting Values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   100         75            7.2462237576E+01  1.7340283401E+00
  b2 =     1          2.5          2.6180768402E+00  8.8295217536E-02
  b3 =     0.1        0.07         6.7359200066E-02  3.4465663377E-03

Residual Sum of Squares:                    8.0565229338E+00
Residual Standard Deviation:                1.1587725499E+00
Degrees of Freedom:                                6
Number of Observations:                            9

$offtext

*-----------------------------------------------------------------------------
* data
*-----------------------------------------------------------------------------

set i /i1*i9/;

table data(i,*)
                  y               x
  i1           8.930E0        9.000E0
  i2          10.800E0       14.000E0
  i3          18.590E0       21.000E0
  i4          22.330E0       28.000E0
  i5          39.350E0       42.000E0
  i6          56.110E0       57.000E0
  i7          61.730E0       63.000E0
  i8          64.620E0       70.000E0
  i9          67.080E0       79.000E0

;

*
* extract data
*
parameter x(i),y(i);
x(i) = data(i,'x');
y(i) = data(i,'y');

*
* certified values
*
scalars
  cb1 'certified value for b1' / 7.2462237576E+01 /
  cb2 'certified value for b2' / 2.6180768402E+00 /
  cb3 'certified value for b3' / 6.7359200066E-02 /
  ce1 'certified std err for b1 ' / 1.7340283401E+00 /
  ce2 'certified std err for b2 ' / 8.8295217536E-02 /
  ce3 'certified std err for b3 ' / 3.4465663377E-03 /
;

*-----------------------------------------------------------------------------
* statistical model
*-----------------------------------------------------------------------------


variables
    sse          'sum of squared errors'
    b1           'coefficient to estimate'
    b2           'coefficient to estimate'
    b3           'coefficient to estimate'
;

equations
    fit(i)     'the non-linear model'
    obj        'objective'
;

obj..     sse =n= 0;
fit(i)..  y(i) =e= b1 / (1+exp[b2-b3*x(i)]);

option nlp=nls;
model nlfit /obj,fit/;

*-----------------------------------------------------------------------------
* first set of initial values
*-----------------------------------------------------------------------------

b1.l =  100;
b2.l =    1;
b3.l =    0.1;

solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.l,b3.l;

abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3))>0.0001) "Accuracy problem";
abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3))>0.0001) "Accuracy problem";

*-----------------------------------------------------------------------------
* second set of initial values
*-----------------------------------------------------------------------------

b1.l =  75;
b2.l =   2.5;
b3.l =   0.07;

solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.l,b3.l;

abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3))>0.0001) "Accuracy problem";
abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3))>0.0001) "Accuracy problem";