$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 circular interference transmittance. The response variable is transmittance, and the predictor variable is wavelength. Reference: Eckerle, K., NIST (197?). Circular Interference Transmittance Study. Data: 1 Response Variable (y = transmittance) 1 Predictor Variable (x = wavelength) 35 Observations Higher Level of Difficulty Observed Data Model: Exponential Class 3 Parameters (b1 to b3) y = (b1/b2) * exp[-0.5*((x-b3)/b2)**2] + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 1 1.5 1.5543827178E+00 1.5408051163E-02 b2 = 10 5 4.0888321754E+00 4.6803020753E-02 b3 = 500 450 4.5154121844E+02 4.6800518816E-02 Residual Sum of Squares: 1.4635887487E-03 Residual Standard Deviation: 6.7629245447E-03 Degrees of Freedom: 32 Number of Observations: 35 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i35/; table data(i,*) y x i1 0.0001575E0 400.000000E0 i2 0.0001699E0 405.000000E0 i3 0.0002350E0 410.000000E0 i4 0.0003102E0 415.000000E0 i5 0.0004917E0 420.000000E0 i6 0.0008710E0 425.000000E0 i7 0.0017418E0 430.000000E0 i8 0.0046400E0 435.000000E0 i9 0.0065895E0 436.500000E0 i10 0.0097302E0 438.000000E0 i11 0.0149002E0 439.500000E0 i12 0.0237310E0 441.000000E0 i13 0.0401683E0 442.500000E0 i14 0.0712559E0 444.000000E0 i15 0.1264458E0 445.500000E0 i16 0.2073413E0 447.000000E0 i17 0.2902366E0 448.500000E0 i18 0.3445623E0 450.000000E0 i19 0.3698049E0 451.500000E0 i20 0.3668534E0 453.000000E0 i21 0.3106727E0 454.500000E0 i22 0.2078154E0 456.000000E0 i23 0.1164354E0 457.500000E0 i24 0.0616764E0 459.000000E0 i25 0.0337200E0 460.500000E0 i26 0.0194023E0 462.000000E0 i27 0.0117831E0 463.500000E0 i28 0.0074357E0 465.000000E0 i29 0.0022732E0 470.000000E0 i30 0.0008800E0 475.000000E0 i31 0.0004579E0 480.000000E0 i32 0.0002345E0 485.000000E0 i33 0.0001586E0 490.000000E0 i34 0.0001143E0 495.000000E0 i35 0.0000710E0 500.000000E0 ; * * 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.5543827178E+00 / cb2 'certified value for b2' / 4.0888321754E+00 / cb3 'certified value for b3' / 4.5154121844E+02 / ce1 'certified std err for b1 ' / 1.5408051163E-02 / ce2 'certified std err for b2 ' / 4.6803020753E-02 / ce3 'certified std err for b3 ' / 4.6800518816E-02 / ; *----------------------------------------------------------------------------- * 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= (b1/b2) * exp[-0.5*sqr((x(i)-b3)/b2)]; option nlp=nls; model nlfit /obj,fit/; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 1; b2.l = 10; b3.l = 500; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.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 = 1.5; b2.l = 5; b3.l = 450; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.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";