$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 the dry weight of onion bulbs
               and tops, and the predictor variable is growing time.

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

Data:          1 Response  (y = onion bulb dry weight)
               1 Predictor (x = growing time)
               15 Observations
               Higher Level of Difficulty
               Observed Data

Model:         Exponential Class
               4 Parameters (b1 to b4)

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

          Starting Values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   100         700           6.9964151270E+02  1.6302297817E+01
  b2 =    10           5           5.2771253025E+00  2.0828735829E+00
  b3 =     1           0.75        7.5962938329E-01  1.9566123451E-01
  b4 =     1           1.3         1.2792483859E+00  6.8761936385E-01

Residual Sum of Squares:                    8.7864049080E+03
Residual Standard Deviation:                2.8262414662E+01
Degrees of Freedom:                                9
Number of Observations:                           15

$offtext

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

set i /i1*i15/;

table data(i,*)
                  y         x
  i1          16.08E0     1.0E0
  i2          33.83E0     2.0E0
  i3          65.80E0     3.0E0
  i4          97.20E0     4.0E0
  i5         191.55E0     5.0E0
  i6         326.20E0     6.0E0
  i7         386.87E0     7.0E0
  i8         520.53E0     8.0E0
  i9         590.03E0     9.0E0
  i10        651.92E0    10.0E0
  i11        724.93E0    11.0E0
  i12        699.56E0    12.0E0
  i13        689.96E0    13.0E0
  i14        637.56E0    14.0E0
  i15        717.41E0    15.0E0

;



*
* 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' / 6.9964151270E+02 /
  cb2 'certified value for b2' / 5.2771253025E+00 /
  cb3 'certified value for b3' / 7.5962938329E-01 /
  cb4 'certified value for b4' / 1.2792483859E+00 /
  ce1 'certified std err for b1 ' / 1.6302297817E+01 /
  ce2 'certified std err for b2 ' / 2.0828735829E+00 /
  ce3 'certified std err for b3 ' / 1.9566123451E-01 /
  ce4 'certified std err for b4 ' / 6.8761936385E-01 /
;

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


variables
    sse          'sum of squared errors'
    b1           'coefficient to estimate'
    b2           'coefficient to estimate'
    b3           'coefficient to estimate'
    b4           '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)])**(1/b4)) ;

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

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

b1.l =  100;
b2.l =   10;
b3.l =    1;
b4.l =    1;

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

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

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

b1.l =   700;
b2.l =     5;
b3.l =     0.75;
b4.l =     1.3;

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

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