$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 the result of a NIST study involving
               ultrasonic calibration.  The response variable is
               ultrasonic response, and the predictor variable is
               metal distance.

Reference:     Chwirut, D., NIST (197?).
               Ultrasonic Reference Block Study.

Data:          1 Response  (y = ultrasonic response)
               1 Predictor (x = metal distance)
               54 Observations
               Lower Level of Difficulty
               Observed Data

Model:         Exponential Class
               3 Parameters (b1 to b3)

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

          Starting values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   0.1         0.15          1.6657666537E-01  3.8303286810E-02
  b2 =   0.01        0.008         5.1653291286E-03  6.6621605126E-04
  b3 =   0.02        0.010         1.2150007096E-02  1.5304234767E-03

Residual Sum of Squares:                    5.1304802941E+02
Residual Standard Deviation:                3.1717133040E+00
Degrees of Freedom:                                51
Number of Observations:                            54

$offtext

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

set i /i1*i54/;

table data(i,*)
              y           x
  i1     92.9000E0     0.500E0
  i2     57.1000E0     1.000E0
  i3     31.0500E0     1.750E0
  i4     11.5875E0     3.750E0
  i5      8.0250E0     5.750E0
  i6     63.6000E0     0.875E0
  i7     21.4000E0     2.250E0
  i8     14.2500E0     3.250E0
  i9      8.4750E0     5.250E0
  i10    63.8000E0     0.750E0
  i11    26.8000E0     1.750E0
  i12    16.4625E0     2.750E0
  i13     7.1250E0     4.750E0
  i14    67.3000E0     0.625E0
  i15    41.0000E0     1.250E0
  i16    21.1500E0     2.250E0
  i17     8.1750E0     4.250E0
  i18    81.5000E0      .500E0
  i19    13.1200E0     3.000E0
  i20    59.9000E0      .750E0
  i21    14.6200E0     3.000E0
  i22    32.9000E0     1.500E0
  i23     5.4400E0     6.000E0
  i24    12.5600E0     3.000E0
  i25     5.4400E0     6.000E0
  i26    32.0000E0     1.500E0
  i27    13.9500E0     3.000E0
  i28    75.8000E0      .500E0
  i29    20.0000E0     2.000E0
  i30    10.4200E0     4.000E0
  i31    59.5000E0      .750E0
  i32    21.6700E0     2.000E0
  i33     8.5500E0     5.000E0
  i34    62.0000E0      .750E0
  i35    20.2000E0     2.250E0
  i36     7.7600E0     3.750E0
  i37     3.7500E0     5.750E0
  i38    11.8100E0     3.000E0
  i39    54.7000E0      .750E0
  i40    23.7000E0     2.500E0
  i41    11.5500E0     4.000E0
  i42    61.3000E0      .750E0
  i43    17.7000E0     2.500E0
  i44     8.7400E0     4.000E0
  i45    59.2000E0      .750E0
  i46    16.3000E0     2.500E0
  i47     8.6200E0     4.000E0
  i48    81.0000E0      .500E0
  i49     4.8700E0     6.000E0
  i50    14.6200E0     3.000E0
  i51    81.7000E0      .500E0
  i52    17.1700E0     2.750E0
  i53    81.3000E0      .500E0
  i54    28.9000E0     1.750E0
;

*
* 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' / 1.6657666537E-01 /
  cb2 'certified value for b2' / 5.1653291286E-03 /
  cb3 'certified value for b3' / 1.2150007096E-02 /
  ce1 'certified std err for b1 ' / 3.8303286810E-02 /
  ce2 'certified std err for b2 ' / 6.6621605126E-04 /
  ce3 'certified std err for b3 ' / 1.5304234767E-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= exp(-b1*x(i))/(b2+b3*x(i));

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

b1.l = 0.1;
b2.l = 0.01;
b3.l = 0.02;

option nlp=nls;
model nlfit /obj,fit/;
solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.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 = 0.15;
b2.l = 0.008;
b3.l = 0.010;

solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.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";