$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 taken from an example discussed in
               Lanczos (1956).  The data were generated to 5-digits
               of accuracy using
               f(x) = 0.0951*exp(-x) + 0.8607*exp(-3*x)
                                     + 1.5576*exp(-5*x).


Reference:     Lanczos, C. (1956).
               Applied Analysis.
               Englewood Cliffs, NJ:  Prentice Hall, pp. 272-280.




Data:          1 Response  (y)
               1 Predictor (x)
               24 Observations
               Lower Level of Difficulty
               Generated Data

Model:         Exponential Class
               6 Parameters (b1 to b6)

               y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x)  +  e



          Starting values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   1.2         0.5           8.6816414977E-02  1.7197908859E-02
  b2 =   0.3         0.7           9.5498101505E-01  9.7041624475E-02
  b3 =   5.6         3.6           8.4400777463E-01  4.1488663282E-02
  b4 =   5.5         4.2           2.9515951832E+00  1.0766312506E-01
  b5 =   6.5         4             1.5825685901E+00  5.8371576281E-02
  b6 =   7.6         6.3           4.9863565084E+00  3.4436403035E-02

Residual Sum of Squares:                    1.6117193594E-08
Residual Standard Deviation:                2.9923229172E-05
Degrees of Freedom:                                18
Number of Observations:                            24

$offtext

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

set i /i1*i24/;

table data(i,*)

               y           x
  i1      2.5134E+00  0.00000E+00
  i2      2.0443E+00  5.00000E-02
  i3      1.6684E+00  1.00000E-01
  i4      1.3664E+00  1.50000E-01
  i5      1.1232E+00  2.00000E-01
  i6      0.9269E+00  2.50000E-01
  i7      0.7679E+00  3.00000E-01
  i8      0.6389E+00  3.50000E-01
  i9      0.5338E+00  4.00000E-01
  i10     0.4479E+00  4.50000E-01
  i11     0.3776E+00  5.00000E-01
  i12     0.3197E+00  5.50000E-01
  i13     0.2720E+00  6.00000E-01
  i14     0.2325E+00  6.50000E-01
  i15     0.1997E+00  7.00000E-01
  i16     0.1723E+00  7.50000E-01
  i17     0.1493E+00  8.00000E-01
  i18     0.1301E+00  8.50000E-01
  i19     0.1138E+00  9.00000E-01
  i20     0.1000E+00  9.50000E-01
  i21     0.0883E+00  1.00000E+00
  i22     0.0783E+00  1.05000E+00
  i23     0.0698E+00  1.10000E+00
  i24     0.0624E+00  1.15000E+00

;

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

set j /j1*j3/;

*
* certified values
*
parameters
  ca(j) 'certified value for a' /j1 8.6816414977E-02
                                 j2 8.4400777463E-01
                                 j3 1.5825685901E+00  /
  cb(j) 'certified value for b' /j1 9.5498101505E-01
                                 j2 2.9515951832E+00
                                 j3 4.9863565084E+00 /
  cea(j) 'certified std err for a ' /j1 1.7197908859E-02
                                     j2 4.1488663282E-02
                                     j3 5.8371576281E-02 /
  ceb(j) 'certified std err for b ' /j1 9.7041624475E-02
                                     j2 1.0766312506E-01
                                     j3 3.4436403035E-02 /
;

set st 'start value' /st1*st3/;

*
* starting values
*
table start(*,j,st)
                       st1               st2               st3
  a.j1    1.2000000000E+00  5.0000000000E-01  8.6816414977E-02
  a.j2    5.6000000000E+00  3.6000000000E+00  8.4400777463E-01
  a.j3    6.5000000000E+00  4.0000000000E+00  1.5825685901E+00
  b.j1    3.0000000000E-01  7.0000000000E-01  9.5498101505E-01
  b.j2    5.5000000000E+00  4.2000000000E+00  2.9515951832E+00
  b.j3    7.6000000000E+00  6.3000000000E+00  4.9863565084E+00
;

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


variables
    sse          'sum of squared errors'
    a(j)         'coefficient to estimate'
    b(j)         'coefficient to estimate'
;

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

obj..     sse =n= 0;
fit(i)..  y(i) =e= sum(j, a(j)*exp[-b(j)*x(i)] );

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

*-----------------------------------------------------------------------------
* solve with different starting points
*-----------------------------------------------------------------------------

loop(st,

   a.l(j) = start("a",j,st);
   b.l(j) = start("b",j,st);

   solve nlfit minimizing sse using nlp;

   abort$(sum(j,abs[a.l(j)-ca(j)]+abs[b.l(j)-cb(j)])>0.0001) "Accuracy problem";
   abort$(sum(j,abs[a.m(j)-cea(j)]+abs[b.m(j)-ceb(j)])>0.0001) "Accuracy problem";

);