$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 semiconductor electron mobility. The response variable is a measure of electron mobility, and the predictor variable is the natural log of the density. Reference: Thurber, R., NIST (197?). Semiconductor electron mobility modeling. Data: 1 Response Variable (y = electron mobility) 1 Predictor Variable (x = log[density]) 37 Observations Higher Level of Difficulty Observed Data Model: Rational Class (cubic/cubic) 7 Parameters (b1 to b7) y = (b1 + b2*x + b3*x**2 + b4*x**3) / (1 + b5*x + b6*x**2 + b7*x**3) + e Starting Values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 1000 1300 1.2881396800E+03 4.6647963344E+00 b2 = 1000 1500 1.4910792535E+03 3.9571156086E+01 b3 = 400 500 5.8323836877E+02 2.8698696102E+01 b4 = 40 75 7.5416644291E+01 5.5675370270E+00 b5 = 0.7 1 9.6629502864E-01 3.1333340687E-02 b6 = 0.3 0.4 3.9797285797E-01 1.4984928198E-02 b7 = 0.03 0.05 4.9727297349E-02 6.5842344623E-03 Residual Sum of Squares: 5.6427082397E+03 Residual Standard Deviation: 1.3714600784E+01 Degrees of Freedom: 30 Number of Observations: 37 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i37/; table data(i,*) y x i1 80.574E0 -3.067E0 i2 84.248E0 -2.981E0 i3 87.264E0 -2.921E0 i4 87.195E0 -2.912E0 i5 89.076E0 -2.840E0 i6 89.608E0 -2.797E0 i7 89.868E0 -2.702E0 i8 90.101E0 -2.699E0 i9 92.405E0 -2.633E0 i10 95.854E0 -2.481E0 i11 100.696E0 -2.363E0 i12 101.060E0 -2.322E0 i13 401.672E0 -1.501E0 i14 390.724E0 -1.460E0 i15 567.534E0 -1.274E0 i16 635.316E0 -1.212E0 i17 733.054E0 -1.100E0 i18 759.087E0 -1.046E0 i19 894.206E0 -0.915E0 i20 990.785E0 -0.714E0 i21 1090.109E0 -0.566E0 i22 1080.914E0 -0.545E0 i23 1122.643E0 -0.400E0 i24 1178.351E0 -0.309E0 i25 1260.531E0 -0.109E0 i26 1273.514E0 -0.103E0 i27 1288.339E0 0.010E0 i28 1327.543E0 0.119E0 i29 1353.863E0 0.377E0 i30 1414.509E0 0.790E0 i31 1425.208E0 0.963E0 i32 1421.384E0 1.006E0 i33 1442.962E0 1.115E0 i34 1464.350E0 1.572E0 i35 1468.705E0 1.841E0 i36 1447.894E0 2.047E0 i37 1457.628E0 2.200E0 ; * * 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.2881396800E+03 / cb2 'certified value for b2' / 1.4910792535E+03 / cb3 'certified value for b3' / 5.8323836877E+02 / cb4 'certified value for b4' / 7.5416644291E+01 / cb5 'certified value for b5' / 9.6629502864E-01 / cb6 'certified value for b6' / 3.9797285797E-01 / cb7 'certified value for b7' / 4.9727297349E-02 / ce1 'certified std err for b1 ' / 4.6647963344E+00 / ce2 'certified std err for b2 ' / 3.9571156086E+01 / ce3 'certified std err for b3 ' / 2.8698696102E+01 / ce4 'certified std err for b4 ' / 5.5675370270E+00 / ce5 'certified std err for b5 ' / 3.1333340687E-02 / ce6 'certified std err for b6 ' / 1.4984928198E-02 / ce7 'certified std err for b7 ' / 6.5842344623E-03 / ; *----------------------------------------------------------------------------- * statistical model *----------------------------------------------------------------------------- variables sse 'sum of squared errors' b1 'coefficient to estimate' b2 'coefficient to estimate' b3 'coefficient to estimate' b4 'coefficient to estimate' b5 'coefficient to estimate' b6 'coefficient to estimate' b7 '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) + b3*sqr(x(i)) + b4*power(x(i),3)) / (1 + b5*x(i) + b6*sqr(x(i)) + b7*power(x(i),3)); option nlp=nls; model nlfit /obj,fit/; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 1000; b2.l = 1000; b3.l = 400; b4.l = 40; b5.l = 0.7; b6.l = 0.3; b7.l = 0.03; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l,b6.l,b7.l; abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3) +abs(b4.l-cb4)+abs(b5.l-cb5)+abs(b6.l-cb6) +abs(b7.l-cb7))>0.0001) "Accuracy problem"; abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3) +abs(b4.m-ce4)+abs(b5.m-ce5)+abs(b6.m-ce6) +abs(b7.m-ce7))>0.0001) "Accuracy problem"; *----------------------------------------------------------------------------- * second set of initial values *----------------------------------------------------------------------------- b1.l = 1300; b2.l = 1500; b3.l = 500; b4.l = 75; b5.l = 1; b6.l = 0.4; b7.l = 0.05; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l,b6.l,b7.l; abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3) +abs(b4.l-cb4)+abs(b5.l-cb5)+abs(b6.l-cb6) +abs(b7.l-cb7))>0.0001) "Accuracy problem"; abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3) +abs(b4.m-ce4)+abs(b5.m-ce5)+abs(b6.m-ce6) +abs(b7.m-ce7))>0.0001) "Accuracy problem";