$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:   These data are described in detail in Box, Hunter and
               Hunter (1978).  The response variable is biochemical
               oxygen demand (BOD) in mg/l, and the predictor
               variable is incubation time in days.

Reference:     Box, G. P., W. G. Hunter, and J. S. Hunter (1978).
               Statistics for Experimenters.
               New York, NY: Wiley, pp. 483-487.

Data:          1 Response  (y = biochemical oxygen demand)
               1 Predictor (x = incubation time)
               6 Observations
               Higher Level of Difficulty
               Observed Data

Model:         Exponential Class
               2 Parameters (b1 and b2)

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

          Starting values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   1           100           2.1380940889E+02  1.2354515176E+01
  b2 =   1             0.75        5.4723748542E-01  1.0455993237E-01

Residual Sum of Squares:                    1.1680088766E+03
Residual Standard Deviation:                1.7088072423E+01
Degrees of Freedom:                                4
Number of Observations:                            6

$offtext

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

set i /i1*i6/;

table data(i,*)
           y     x
  i1     109     1
  i2     149     2
  i3     149     3
  i4     191     5
  i5     213     7
  i6     224    10
;


*
* 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' / 2.1380940889E+02 /
  cb2 'certified value for b2' / 5.4723748542E-01 /
  ce1 'certified std err for b1 ' / 1.2354515176E+01 /
  ce2 'certified std err for b2 ' / 1.0455993237E-01 /
;


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


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

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

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

$ontext

does not converge

b1.l = 1;
b2.l = 1;

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

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

$offtext
*-----------------------------------------------------------------------------
* second set of initial values
*-----------------------------------------------------------------------------

b1.l =  100;
b2.l =    0.75;

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

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