$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
               quantum defects in iodine atoms.  The response
               variable is the number of quantum defects, and the
               predictor variable is the excited energy state.
               The argument to the ARCTAN function is in radians.

Reference:     Roszman, L., NIST (19??).
               Quantum Defects for Sulfur I Atom.

Data:          1 Response  (y = quantum defect)
               1 Predictor (x = excited state energy)
               25 Observations
               Average Level of Difficulty
               Observed Data

Model:         Miscellaneous Class
               4 Parameters (b1 to b4)

               pi = 3.141592653589793238462643383279E0
               y =  b1 - b2*x - arctan[b3/(x-b4)]/pi  +  e


          Starting Values                  Certified Values

        Start 1     Start 2           Parameter     Standard Deviation
  b1 =      0.1         0.2         2.0196866396E-01  1.9172666023E-02
  b2 =     -0.00001    -0.000005   -6.1953516256E-06  3.2058931691E-06
  b3 =   1000        1200           1.2044556708E+03  7.4050983057E+01
  b4 =   -100        -150          -1.8134269537E+02  4.9573513849E+01

Residual Sum of Squares:                    4.9484847331E-04
Residual Standard Deviation:                4.8542984060E-03
Degrees of Freedom:                                 21
Number of Observations:                             25
$offtext

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

set i /i1*i25/;

table data(i,*)
             y           x
  i1     0.252429    -4868.68
  i2     0.252141    -4868.09
  i3     0.251809    -4867.41
  i4     0.297989    -3375.19
  i5     0.296257    -3373.14
  i6     0.295319    -3372.03
  i7     0.339603    -2473.74
  i8     0.337731    -2472.35
  i9     0.333820    -2469.45
  i10    0.389510    -1894.65
  i11    0.386998    -1893.40
  i12    0.438864    -1497.24
  i13    0.434887    -1495.85
  i14    0.427893    -1493.41
  i15    0.471568    -1208.68
  i16    0.461699    -1206.18
  i17    0.461144    -1206.04
  i18    0.513532     -997.92
  i19    0.506641     -996.61
  i20    0.505062     -996.31
  i21    0.535648     -834.94
  i22    0.533726     -834.66
  i23    0.568064     -710.03
  i24    0.612886     -530.16
  i25    0.624169     -464.17

;

*
* 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' /  2.0196866396E-01 /
  cb2 'certified value for b2' / -6.1953516256E-06 /
  cb3 'certified value for b3' /  1.2044556708E+03 /
  cb4 'certified value for b4' / -1.8134269537E+02 /
  ce1 'certified std err for b1 ' / 1.9172666023E-02 /
  ce2 'certified std err for b2 ' / 3.2058931691E-06 /
  ce3 'certified std err for b3 ' / 7.4050983057E+01 /
  ce4 'certified std err for b4 ' / 4.9573513849E+01 /
;



*-----------------------------------------------------------------------------
* 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 - b2*x(i) - arctan[b3/(x(i)-b4)]/pi;


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

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

b1.l =     0.1;
b2.l =    -0.00001;
b3.l =  1000;
b4.l =  -100;

solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.l,b3.l,b4.l;

abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3)
       +abs(b4.l-cb4))>0.001) "Accuracy problem";
abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3)
       +abs(b4.m-ce4))>0.001) "Accuracy problem";


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


b1.l =     0.2;
b2.l =    -0.000005;
b3.l =  1200;
b4.l =  -150;

solve nlfit minimizing sse using nlp;
display sse.l,b1.l,b2.l,b3.l,b4.l;

abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3)
       +abs(b4.l-cb4))>0.001) "Accuracy problem";
abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3)
       +abs(b4.m-ce4))>0.001) "Accuracy problem";