$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 and model are described in Daniel and Wood
               (1980), and originally published in E.S.Keeping,
               "Introduction to Statistical Inference," Van Nostrand
               Company, Princeton, NJ, 1962, p. 354.  The response
               variable is energy radieted from a carbon filament
               lamp per cm**2 per second, and the predictor variable
               is the absolute temperature of the filament in 1000
               degrees Kelvin.

Reference:     Daniel, C. and F. S. Wood (1980).
               Fitting Equations to Data, Second Edition.
               New York, NY:  John Wiley and Sons, pp. 428-431.

Data:          1 Response Variable  (y = energy)
               1 Predictor Variable (x = temperature)
               6 Observations
               Lower Level of Difficulty
               Observed Data

Model:         Miscellaneous Class
               2 Parameters (b1 and b2)

               y  = b1*x**b2  +  e

          Starting values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   1           0.7           7.6886226176E-01  1.8281973860E-02
  b2 =   5           4             3.8604055871E+00  5.1726610913E-02

Residual Sum of Squares:                    4.3173084083E-03
Residual Standard Deviation:                3.2853114039E-02
Degrees of Freedom:                                4
Number of Observations:                            6

$offtext


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

set i /i1*i6/;

table data(i,*)
             y           x
  i1   2.138E0     1.309E0
  i2   3.421E0     1.471E0
  i3   3.597E0     1.490E0
  i4   4.340E0     1.565E0
  i5   4.882E0     1.611E0
  i6   5.660E0     1.680E0
;

*
* 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' /7.6886226176E-01/
  cb2 'certified value for b2' /3.8604055871E+00/
  ce1 'certified std err for b1 ' / 1.8281973860E-02 /
  ce2 'certified std err for b2 ' / 5.1726610913E-02 /
;

*-----------------------------------------------------------------------------
* 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*x(i)**b2;

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

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

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))>0.0001) "Accuracy problem";
abort$((abs(b1.m-ce1)+abs(b2.m-ce2))>0.0001) "Accuracy problem";

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

b1.l = 0.7;
b2.l = 4;

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";

*-----------------------------------------------------------------------------
* third set of initial values
*-----------------------------------------------------------------------------

b1.l = 7.6886226176E-01;
b2.l = 3.8604055871E+00;

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";