$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 superconductivity magnetization modeling. The response variable is magnetism, and the predictor variable is the log of time in minutes. Reference: Bennett, L., L. Swartzendruber, and H. Brown, NIST (1994). Superconductivity Magnetization Modeling. Data: 1 Response Variable (y = magnetism) 1 Predictor Variable (x = log[time]) 154 Observations Higher Level of Difficulty Observed Data Model: Miscellaneous Class 3 Parameters (b1 to b3) y = b1 * (b2+x)**(-1/b3) + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = -2000 -1500 -2.5235058043E+03 2.9715175411E+02 b2 = 50 45 4.6736564644E+01 1.2448871856E+00 b3 = 0.8 0.85 9.3218483193E-01 2.0272299378E-02 Residual Sum of Squares: 5.2404744073E-04 Residual Standard Deviation: 1.8629312528E-03 Degrees of Freedom: 151 Number of Observations: 154 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i154/; table data(i,*) y x i1 -34.834702E0 7.447168E0 i2 -34.393200E0 8.102586E0 i3 -34.152901E0 8.452547E0 i4 -33.979099E0 8.711278E0 i5 -33.845901E0 8.916774E0 i6 -33.732899E0 9.087155E0 i7 -33.640301E0 9.232590E0 i8 -33.559200E0 9.359535E0 i9 -33.486801E0 9.472166E0 i10 -33.423100E0 9.573384E0 i11 -33.365101E0 9.665293E0 i12 -33.313000E0 9.749461E0 i13 -33.260899E0 9.827092E0 i14 -33.217400E0 9.899128E0 i15 -33.176899E0 9.966321E0 i16 -33.139198E0 10.029280E0 i17 -33.101601E0 10.088510E0 i18 -33.066799E0 10.144430E0 i19 -33.035000E0 10.197380E0 i20 -33.003101E0 10.247670E0 i21 -32.971298E0 10.295560E0 i22 -32.942299E0 10.341250E0 i23 -32.916302E0 10.384950E0 i24 -32.890202E0 10.426820E0 i25 -32.864101E0 10.467000E0 i26 -32.841000E0 10.505640E0 i27 -32.817799E0 10.542830E0 i28 -32.797501E0 10.578690E0 i29 -32.774300E0 10.613310E0 i30 -32.757000E0 10.646780E0 i31 -32.733799E0 10.679150E0 i32 -32.716400E0 10.710520E0 i33 -32.699100E0 10.740920E0 i34 -32.678799E0 10.770440E0 i35 -32.661400E0 10.799100E0 i36 -32.644001E0 10.826970E0 i37 -32.626701E0 10.854080E0 i38 -32.612202E0 10.880470E0 i39 -32.597698E0 10.906190E0 i40 -32.583199E0 10.931260E0 i41 -32.568699E0 10.955720E0 i42 -32.554298E0 10.979590E0 i43 -32.539799E0 11.002910E0 i44 -32.525299E0 11.025700E0 i45 -32.510799E0 11.047980E0 i46 -32.499199E0 11.069770E0 i47 -32.487598E0 11.091100E0 i48 -32.473202E0 11.111980E0 i49 -32.461601E0 11.132440E0 i50 -32.435501E0 11.152480E0 i51 -32.435501E0 11.172130E0 i52 -32.426800E0 11.191410E0 i53 -32.412300E0 11.210310E0 i54 -32.400799E0 11.228870E0 i55 -32.392101E0 11.247090E0 i56 -32.380501E0 11.264980E0 i57 -32.366001E0 11.282560E0 i58 -32.357300E0 11.299840E0 i59 -32.348598E0 11.316820E0 i60 -32.339901E0 11.333520E0 i61 -32.328400E0 11.349940E0 i62 -32.319698E0 11.366100E0 i63 -32.311001E0 11.382000E0 i64 -32.299400E0 11.397660E0 i65 -32.290699E0 11.413070E0 i66 -32.282001E0 11.428240E0 i67 -32.273300E0 11.443200E0 i68 -32.264599E0 11.457930E0 i69 -32.256001E0 11.472440E0 i70 -32.247299E0 11.486750E0 i71 -32.238602E0 11.500860E0 i72 -32.229900E0 11.514770E0 i73 -32.224098E0 11.528490E0 i74 -32.215401E0 11.542020E0 i75 -32.203800E0 11.555380E0 i76 -32.198002E0 11.568550E0 i77 -32.189400E0 11.581560E0 i78 -32.183601E0 11.594420E0 i79 -32.174900E0 11.607121E0 i80 -32.169102E0 11.619640E0 i81 -32.163300E0 11.632000E0 i82 -32.154598E0 11.644210E0 i83 -32.145901E0 11.656280E0 i84 -32.140099E0 11.668200E0 i85 -32.131401E0 11.679980E0 i86 -32.125599E0 11.691620E0 i87 -32.119801E0 11.703130E0 i88 -32.111198E0 11.714510E0 i89 -32.105400E0 11.725760E0 i90 -32.096699E0 11.736880E0 i91 -32.090900E0 11.747890E0 i92 -32.088001E0 11.758780E0 i93 -32.079300E0 11.769550E0 i94 -32.073502E0 11.780200E0 i95 -32.067699E0 11.790730E0 i96 -32.061901E0 11.801160E0 i97 -32.056099E0 11.811480E0 i98 -32.050301E0 11.821700E0 i99 -32.044498E0 11.831810E0 i100 -32.038799E0 11.841820E0 i101 -32.033001E0 11.851730E0 i102 -32.027199E0 11.861550E0 i103 -32.024300E0 11.871270E0 i104 -32.018501E0 11.880890E0 i105 -32.012699E0 11.890420E0 i106 -32.004002E0 11.899870E0 i107 -32.001099E0 11.909220E0 i108 -31.995300E0 11.918490E0 i109 -31.989500E0 11.927680E0 i110 -31.983700E0 11.936780E0 i111 -31.977900E0 11.945790E0 i112 -31.972099E0 11.954730E0 i113 -31.969299E0 11.963590E0 i114 -31.963501E0 11.972370E0 i115 -31.957701E0 11.981070E0 i116 -31.951900E0 11.989700E0 i117 -31.946100E0 11.998260E0 i118 -31.940300E0 12.006740E0 i119 -31.937401E0 12.015150E0 i120 -31.931601E0 12.023490E0 i121 -31.925800E0 12.031760E0 i122 -31.922899E0 12.039970E0 i123 -31.917101E0 12.048100E0 i124 -31.911301E0 12.056170E0 i125 -31.908400E0 12.064180E0 i126 -31.902599E0 12.072120E0 i127 -31.896900E0 12.080010E0 i128 -31.893999E0 12.087820E0 i129 -31.888201E0 12.095580E0 i130 -31.885300E0 12.103280E0 i131 -31.882401E0 12.110920E0 i132 -31.876600E0 12.118500E0 i133 -31.873699E0 12.126030E0 i134 -31.867901E0 12.133500E0 i135 -31.862101E0 12.140910E0 i136 -31.859200E0 12.148270E0 i137 -31.856300E0 12.155570E0 i138 -31.850500E0 12.162830E0 i139 -31.844700E0 12.170030E0 i140 -31.841801E0 12.177170E0 i141 -31.838900E0 12.184270E0 i142 -31.833099E0 12.191320E0 i143 -31.830200E0 12.198320E0 i144 -31.827299E0 12.205270E0 i145 -31.821600E0 12.212170E0 i146 -31.818701E0 12.219030E0 i147 -31.812901E0 12.225840E0 i148 -31.809999E0 12.232600E0 i149 -31.807100E0 12.239320E0 i150 -31.801300E0 12.245990E0 i151 -31.798401E0 12.252620E0 i152 -31.795500E0 12.259200E0 i153 -31.789700E0 12.265750E0 i154 -31.786800E0 12.272240E0 ; * * 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.5235058043E+03 / cb2 'certified value for b2' / 4.6736564644E+01 / cb3 'certified value for b3' / 9.3218483193E-01 / ce1 'certified std err for b1 ' / 2.9715175411E+02 / ce2 'certified std err for b2 ' / 1.2448871856E+00 / ce3 'certified std err for b3 ' / 2.0272299378E-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+x(i))**(-1/b3) ; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- * We need to extend the max number of function evaluation calls. $onecho > nls.opt * This option is needed in Bennett5.gms mxfcal 1000 $offecho b1.l = -2000; b2.l = 50; b3.l = 0.8; option nlp=nls; model nlfit /obj,fit/; nlfit.optfile=1; 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))>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 = -1500; b2.l = 45; b3.l = 0.85; 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))>0.0001) "Accuracy problem"; abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3))>0.0001) "Accuracy problem";