$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 14-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.5100000027E-02 5.3347304234E-11 b2 = 0.3 0.7 1.0000000001E+00 2.7473038179E-10 b3 = 5.6 3.6 8.6070000013E-01 1.3576062225E-10 b4 = 5.5 4.2 3.0000000002E+00 3.3308253069E-10 b5 = 6.5 4 1.5575999998E+00 1.8815731448E-10 b6 = 7.6 6.3 5.0000000001E+00 1.1057500538E-10 Residual Sum of Squares: 1.4307867721E-25 Residual Standard Deviation: 8.9156129349E-14 Degrees of Freedom: 18 Number of Observations: 24 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i24/; table data(i,*) y x i1 2.513400000000E+00 0.000000000000E+00 i2 2.044333373291E+00 5.000000000000E-02 i3 1.668404436564E+00 1.000000000000E-01 i4 1.366418021208E+00 1.500000000000E-01 i5 1.123232487372E+00 2.000000000000E-01 i6 9.268897180037E-01 2.500000000000E-01 i7 7.679338563728E-01 3.000000000000E-01 i8 6.388775523106E-01 3.500000000000E-01 i9 5.337835317402E-01 4.000000000000E-01 i10 4.479363617347E-01 4.500000000000E-01 i11 3.775847884350E-01 5.000000000000E-01 i12 3.197393199326E-01 5.500000000000E-01 i13 2.720130773746E-01 6.000000000000E-01 i14 2.324965529032E-01 6.500000000000E-01 i15 1.996589546065E-01 7.000000000000E-01 i16 1.722704126914E-01 7.500000000000E-01 i17 1.493405660168E-01 8.000000000000E-01 i18 1.300700206922E-01 8.500000000000E-01 i19 1.138119324644E-01 9.000000000000E-01 i20 1.000415587559E-01 9.500000000000E-01 i21 8.833209084540E-02 1.000000000000E+00 i22 7.833544019350E-02 1.050000000000E+00 i23 6.976693743449E-02 1.100000000000E+00 i24 6.239312536719E-02 1.150000000000E+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.5100000027E-02 j2 8.6070000013E-01 j3 1.5575999998E+00 / cb(j) 'certified value for b' /j1 1.0000000001E+00 j2 3.0000000002E+00 j3 5.0000000001E+00 / cea(j) 'certified std err for a ' /j1 5.3347304234E-11 j2 1.3576062225E-10 j3 1.8815731448E-10 / ceb(j) 'certified std err for b ' /j1 2.7473038179E-10 j2 3.3308253069E-10 j3 1.1057500538E-10 / ; set st 'start value' /st1*st3/; * * starting values * table start(*,j,st) st1 st2 st3 a.j1 1.2000000000E+00 5.0000000000E-01 9.5100000027E-02 a.j2 5.6000000000E+00 3.6000000000E+00 8.6070000013E-01 a.j3 6.5000000000E+00 4.0000000000E+00 1.5575999998E+00 b.j1 3.0000000000E-01 7.0000000000E-01 1.0000000001E+00 b.j2 5.5000000000E+00 4.2000000000E+00 3.0000000002E+00 b.j3 7.6000000000E+00 6.3000000000E+00 5.0000000001E+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"; );