$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 6-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
               Average 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           9.6251029939E-02  6.6770575477E-04
  b2 =   0.3         0.7           1.0057332849E+00  3.3989646176E-03
  b3 =   5.6         3.6           8.6424689056E-01  1.7185846685E-03
  b4 =   5.5         4.2           3.0078283915E+00  4.1707005856E-03
  b5 =   6.5         4             1.5529016879E+00  2.3744381417E-03
  b6 =   7.6         6.3           5.0028798100E+00  1.3958787284E-03

Residual Sum of Squares:                    2.2299428125E-11
Residual Standard Deviation:                1.1130395851E-06
Degrees of Freedom:                                18
Number of Observations:                            24

$offtext

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

set i /i1*i24/;

table data(i,*)

               y            x
  i1    2.51340E+00  0.00000E+00
  i2    2.04433E+00  5.00000E-02
  i3    1.66840E+00  1.00000E-01
  i4    1.36642E+00  1.50000E-01
  i5    1.12323E+00  2.00000E-01
  i6    9.26890E-01  2.50000E-01
  i7    7.67934E-01  3.00000E-01
  i8    6.38878E-01  3.50000E-01
  i9    5.33784E-01  4.00000E-01
  i10   4.47936E-01  4.50000E-01
  i11   3.77585E-01  5.00000E-01
  i12   3.19739E-01  5.50000E-01
  i13   2.72013E-01  6.00000E-01
  i14   2.32497E-01  6.50000E-01
  i15   1.99659E-01  7.00000E-01
  i16   1.72270E-01  7.50000E-01
  i17   1.49341E-01  8.00000E-01
  i18   1.30070E-01  8.50000E-01
  i19   1.13812E-01  9.00000E-01
  i20   1.00042E-01  9.50000E-01
  i21   8.83321E-02  1.00000E+00
  i22   7.83354E-02  1.05000E+00
  i23   6.97669E-02  1.10000E+00
  i24   6.23931E-02  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 9.6251029939E-02
                                 j2 8.6424689056E-01
                                 j3 1.5529016879E+00  /
  cb(j) 'certified value for b' /j1 1.0057332849E+00
                                 j2 3.0078283915E+00
                                 j3 5.0028798100E+00 /
  cea(j) 'certified std err for a ' /j1 6.6770575477E-04
                                     j2 1.7185846685E-03
                                     j3 2.3744381417E-03 /
  ceb(j) 'certified std err for b ' /j1 3.3989646176E-03
                                     j2 4.1707005856E-03
                                     j3 1.3958787284E-03 /
;

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

*
* starting values
*
table start(*,j,st)
                  st1                 st2                st3
  a.j1       1.2000000000E+00   5.0000000000E-01    9.6251029939E-02
  a.j2       5.6000000000E+00   3.6000000000E+00    8.6424689056E-01
  a.j3       6.5000000000E+00   4.0000000000E+00    1.5529016879E+00
  b.j1       3.0000000000E-01   7.0000000000E-01    1.0057332849E+00
  b.j2       5.5000000000E+00   4.2000000000E+00    3.0078283915E+00
  b.j3       7.6000000000E+00   6.3000000000E+00    5.0028798100E+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";

);