$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 the thermal expansion of copper. The response variable is the coefficient of thermal expansion, and the predictor variable is temperature in degrees kelvin. Reference: Hahn, T., NIST (197?). Copper Thermal Expansion Study. Data: 1 Response (y = coefficient of thermal expansion) 1 Predictor (x = temperature, degrees kelvin) 236 Observations Average 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 = 10 1 1.0776351733E+00 1.7070154742E-01 b2 = -1 -0.1 -1.2269296921E-01 1.2000289189E-02 b3 = 0.05 0.005 4.0863750610E-03 2.2508314937E-04 b4 = -0.00001 -0.000001 -1.4262662514E-06 2.7578037666E-07 b5 = -0.05 -0.005 -5.7609940901E-03 2.4712888219E-04 b6 = 0.001 0.0001 2.4053735503E-04 1.0449373768E-05 b7 = -0.000001 -0.0000001 -1.2314450199E-07 1.3027335327E-08 Residual Sum of Squares: 1.5324382854E+00 Residual Standard Deviation: 8.1803852243E-02 Degrees of Freedom: 229 Number of Observations: 236 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i236/; table data(i,*) y x i1 .591E0 24.41E0 i2 1.547E0 34.82E0 i3 2.902E0 44.09E0 i4 2.894E0 45.07E0 i5 4.703E0 54.98E0 i6 6.307E0 65.51E0 i7 7.03E0 70.53E0 i8 7.898E0 75.70E0 i9 9.470E0 89.57E0 i10 9.484E0 91.14E0 i11 10.072E0 96.40E0 i12 10.163E0 97.19E0 i13 11.615E0 114.26E0 i14 12.005E0 120.25E0 i15 12.478E0 127.08E0 i16 12.982E0 133.55E0 i17 12.970E0 133.61E0 i18 13.926E0 158.67E0 i19 14.452E0 172.74E0 i20 14.404E0 171.31E0 i21 15.190E0 202.14E0 i22 15.550E0 220.55E0 i23 15.528E0 221.05E0 i24 15.499E0 221.39E0 i25 16.131E0 250.99E0 i26 16.438E0 268.99E0 i27 16.387E0 271.80E0 i28 16.549E0 271.97E0 i29 16.872E0 321.31E0 i30 16.830E0 321.69E0 i31 16.926E0 330.14E0 i32 16.907E0 333.03E0 i33 16.966E0 333.47E0 i34 17.060E0 340.77E0 i35 17.122E0 345.65E0 i36 17.311E0 373.11E0 i37 17.355E0 373.79E0 i38 17.668E0 411.82E0 i39 17.767E0 419.51E0 i40 17.803E0 421.59E0 i41 17.765E0 422.02E0 i42 17.768E0 422.47E0 i43 17.736E0 422.61E0 i44 17.858E0 441.75E0 i45 17.877E0 447.41E0 i46 17.912E0 448.7E0 i47 18.046E0 472.89E0 i48 18.085E0 476.69E0 i49 18.291E0 522.47E0 i50 18.357E0 522.62E0 i51 18.426E0 524.43E0 i52 18.584E0 546.75E0 i53 18.610E0 549.53E0 i54 18.870E0 575.29E0 i55 18.795E0 576.00E0 i56 19.111E0 625.55E0 i57 .367E0 20.15E0 i58 .796E0 28.78E0 i59 0.892E0 29.57E0 i60 1.903E0 37.41E0 i61 2.150E0 39.12E0 i62 3.697E0 50.24E0 i63 5.870E0 61.38E0 i64 6.421E0 66.25E0 i65 7.422E0 73.42E0 i66 9.944E0 95.52E0 i67 11.023E0 107.32E0 i68 11.87E0 122.04E0 i69 12.786E0 134.03E0 i70 14.067E0 163.19E0 i71 13.974E0 163.48E0 i72 14.462E0 175.70E0 i73 14.464E0 179.86E0 i74 15.381E0 211.27E0 i75 15.483E0 217.78E0 i76 15.59E0 219.14E0 i77 16.075E0 262.52E0 i78 16.347E0 268.01E0 i79 16.181E0 268.62E0 i80 16.915E0 336.25E0 i81 17.003E0 337.23E0 i82 16.978E0 339.33E0 i83 17.756E0 427.38E0 i84 17.808E0 428.58E0 i85 17.868E0 432.68E0 i86 18.481E0 528.99E0 i87 18.486E0 531.08E0 i88 19.090E0 628.34E0 i89 16.062E0 253.24E0 i90 16.337E0 273.13E0 i91 16.345E0 273.66E0 i92 16.388E0 282.10E0 i93 17.159E0 346.62E0 i94 17.116E0 347.19E0 i95 17.164E0 348.78E0 i96 17.123E0 351.18E0 i97 17.979E0 450.10E0 i98 17.974E0 450.35E0 i99 18.007E0 451.92E0 i100 17.993E0 455.56E0 i101 18.523E0 552.22E0 i102 18.669E0 553.56E0 i103 18.617E0 555.74E0 i104 19.371E0 652.59E0 i105 19.330E0 656.20E0 i106 0.080E0 14.13E0 i107 0.248E0 20.41E0 i108 1.089E0 31.30E0 i109 1.418E0 33.84E0 i110 2.278E0 39.70E0 i111 3.624E0 48.83E0 i112 4.574E0 54.50E0 i113 5.556E0 60.41E0 i114 7.267E0 72.77E0 i115 7.695E0 75.25E0 i116 9.136E0 86.84E0 i117 9.959E0 94.88E0 i118 9.957E0 96.40E0 i119 11.600E0 117.37E0 i120 13.138E0 139.08E0 i121 13.564E0 147.73E0 i122 13.871E0 158.63E0 i123 13.994E0 161.84E0 i124 14.947E0 192.11E0 i125 15.473E0 206.76E0 i126 15.379E0 209.07E0 i127 15.455E0 213.32E0 i128 15.908E0 226.44E0 i129 16.114E0 237.12E0 i130 17.071E0 330.90E0 i131 17.135E0 358.72E0 i132 17.282E0 370.77E0 i133 17.368E0 372.72E0 i134 17.483E0 396.24E0 i135 17.764E0 416.59E0 i136 18.185E0 484.02E0 i137 18.271E0 495.47E0 i138 18.236E0 514.78E0 i139 18.237E0 515.65E0 i140 18.523E0 519.47E0 i141 18.627E0 544.47E0 i142 18.665E0 560.11E0 i143 19.086E0 620.77E0 i144 0.214E0 18.97E0 i145 0.943E0 28.93E0 i146 1.429E0 33.91E0 i147 2.241E0 40.03E0 i148 2.951E0 44.66E0 i149 3.782E0 49.87E0 i150 4.757E0 55.16E0 i151 5.602E0 60.90E0 i152 7.169E0 72.08E0 i153 8.920E0 85.15E0 i154 10.055E0 97.06E0 i155 12.035E0 119.63E0 i156 12.861E0 133.27E0 i157 13.436E0 143.84E0 i158 14.167E0 161.91E0 i159 14.755E0 180.67E0 i160 15.168E0 198.44E0 i161 15.651E0 226.86E0 i162 15.746E0 229.65E0 i163 16.216E0 258.27E0 i164 16.445E0 273.77E0 i165 16.965E0 339.15E0 i166 17.121E0 350.13E0 i167 17.206E0 362.75E0 i168 17.250E0 371.03E0 i169 17.339E0 393.32E0 i170 17.793E0 448.53E0 i171 18.123E0 473.78E0 i172 18.49E0 511.12E0 i173 18.566E0 524.70E0 i174 18.645E0 548.75E0 i175 18.706E0 551.64E0 i176 18.924E0 574.02E0 i177 19.1E0 623.86E0 i178 0.375E0 21.46E0 i179 0.471E0 24.33E0 i180 1.504E0 33.43E0 i181 2.204E0 39.22E0 i182 2.813E0 44.18E0 i183 4.765E0 55.02E0 i184 9.835E0 94.33E0 i185 10.040E0 96.44E0 i186 11.946E0 118.82E0 i187 12.596E0 128.48E0 i188 13.303E0 141.94E0 i189 13.922E0 156.92E0 i190 14.440E0 171.65E0 i191 14.951E0 190.00E0 i192 15.627E0 223.26E0 i193 15.639E0 223.88E0 i194 15.814E0 231.50E0 i195 16.315E0 265.05E0 i196 16.334E0 269.44E0 i197 16.430E0 271.78E0 i198 16.423E0 273.46E0 i199 17.024E0 334.61E0 i200 17.009E0 339.79E0 i201 17.165E0 349.52E0 i202 17.134E0 358.18E0 i203 17.349E0 377.98E0 i204 17.576E0 394.77E0 i205 17.848E0 429.66E0 i206 18.090E0 468.22E0 i207 18.276E0 487.27E0 i208 18.404E0 519.54E0 i209 18.519E0 523.03E0 i210 19.133E0 612.99E0 i211 19.074E0 638.59E0 i212 19.239E0 641.36E0 i213 19.280E0 622.05E0 i214 19.101E0 631.50E0 i215 19.398E0 663.97E0 i216 19.252E0 646.9E0 i217 19.89E0 748.29E0 i218 20.007E0 749.21E0 i219 19.929E0 750.14E0 i220 19.268E0 647.04E0 i221 19.324E0 646.89E0 i222 20.049E0 746.9E0 i223 20.107E0 748.43E0 i224 20.062E0 747.35E0 i225 20.065E0 749.27E0 i226 19.286E0 647.61E0 i227 19.972E0 747.78E0 i228 20.088E0 750.51E0 i229 20.743E0 851.37E0 i230 20.83E0 845.97E0 i231 20.935E0 847.54E0 i232 21.035E0 849.93E0 i233 20.93E0 851.61E0 i234 21.074E0 849.75E0 i235 21.085E0 850.98E0 i236 20.935E0 848.23E0 ; * * 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.0776351733E+00 / cb2 'certified value for b2' / -1.2269296921E-01 / cb3 'certified value for b3' / 4.0863750610E-03 / cb4 'certified value for b4' / -1.4262662514E-06 / cb5 'certified value for b5' / -5.7609940901E-03 / cb6 'certified value for b6' / 2.4053735503E-04 / cb7 'certified value for b7' / -1.2314450199E-07 / ce1 'certified std err for b1 ' / 1.7070154742E-01 / ce2 'certified std err for b2 ' / 1.2000289189E-02 / ce3 'certified std err for b3 ' / 2.2508314937E-04 / ce4 'certified std err for b4 ' / 2.7578037666E-07 / ce5 'certified std err for b5 ' / 2.4712888219E-04 / ce6 'certified std err for b6 ' / 1.0449373768E-05 / ce7 'certified std err for b7 ' / 1.3027335327E-08 / ; *----------------------------------------------------------------------------- * 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*x(i)**2 + b4*x(i)**3) / (1 + b5*x(i) + b6*x(i)**2 + b7*x(i)**3); *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 10; b2.l = -1; b3.l = 0.05; b4.l = -0.00001; b5.l = -0.05; b6.l = 0.001; b7.l = -0.000001; option nlp=nls; model nlfit /obj,fit/; 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 = 1; b2.l = -0.1; b3.l = 0.005; b4.l = -0.000001; b5.l = -0.005; b6.l = 0.0001; b7.l = -0.0000001; 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";