$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:   This problem was found to be difficult for some very
               good algorithms.  There is a local minimum at (+inf,
               -14.07..., -inf, -inf) with final sum of squares
               0.00102734....

               See More, J. J., Garbow, B. S., and Hillstrom, K. E.
               (1981).  Testing unconstrained optimization software.
               ACM Transactions on Mathematical Software. 7(1):
               pp. 17-41.

Reference:     Kowalik, J.S., and M. R. Osborne, (1978).
               Methods for Unconstrained Optimization Problems.
               New York, NY:  Elsevier North-Holland.

Data:          1 Response  (y)
               1 Predictor (x)
               11 Observations
               Higher Level of Difficulty
               Generated Data

Model:         Rational Class (linear/quadratic)
               4 Parameters (b1 to b4)

               y = b1*(x**2+x*b2) / (x**2+x*b3+b4)  +  e

          Starting values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =   25          0.25          1.9280693458E-01  1.1435312227E-02
  b2 =   39          0.39          1.9128232873E-01  1.9633220911E-01
  b3 =   41.5        0.415         1.2305650693E-01  8.0842031232E-02
  b4 =   39          0.39          1.3606233068E-01  9.0025542308E-02

Residual Sum of Squares:                    3.0750560385E-04
Residual Standard Deviation:                6.6279236551E-03
Degrees of Freedom:                                7
Number of Observations:                           11

$offtext


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

set i /i1*i11/;

table data(i,*)
             y               x
  i1    1.957000E-01    4.000000E+00
  i2    1.947000E-01    2.000000E+00
  i3    1.735000E-01    1.000000E+00
  i4    1.600000E-01    5.000000E-01
  i5    8.440000E-02    2.500000E-01
  i6    6.270000E-02    1.670000E-01
  i7    4.560000E-02    1.250000E-01
  i8    3.420000E-02    1.000000E-01
  i9    3.230000E-02    8.330000E-02
  i10   2.350000E-02    7.140000E-02
  i11   2.460000E-02    6.250000E-02
;

*
* 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.9280693458E-01/
  cb2 'certified value for b2' /1.9128232873E-01/
  cb3 'certified value for b3' /1.2305650693E-01/
  cb4 'certified value for b4' /1.3606233068E-01/
  ce1 'certified std err for b1 ' / 1.1435312227E-02 /
  ce2 'certified std err for b2 ' / 1.9633220911E-01 /
  ce3 'certified std err for b3 ' / 8.0842031232E-02 /
  ce4 'certified std err for b4 ' / 9.0025542308E-02 /
;

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


variables
    sse          'sum of squared errors'
    b1           'coefficient to estimate'
    b2           'coefficient to estimate'
    b3           'coefficient to estimate'
    b4           'coefficient to estimate'
;

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

obj..     sse =n= 0;
fit(i)..  y(i) =e= b1*(sqr(x(i))+x(i)*b2) / (sqr(x(i))+x(i)*b3+b4);

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

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

b1.l =  25;
b2.l =  39;
b3.l =  41.5;
b4.l =  39;

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

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

b1.l =  0.25;
b2.l =  0.39;
b3.l =  0.415;
b4.l =  0.39;


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