$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: The data are two well-separated Gaussians on a decaying exponential baseline plus normally distributed zero-mean noise with variance = 6.25. Reference: Rust, B., NIST (1996). Data: 1 Response (y) 1 Predictor (x) 250 Observations Lower Level of Difficulty Generated Data Model: Exponential Class 8 Parameters (b1 to b8) y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 ) + b6*exp( -(x-b7)**2 / b8**2 ) + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 97.0 94.0 9.8778210871E+01 5.7527312730E-01 b2 = 0.009 0.0105 1.0497276517E-02 1.1406289017E-04 b3 = 100.0 99.0 1.0048990633E+02 5.8831775752E-01 b4 = 65.0 63.0 6.7481111276E+01 1.0460593412E-01 b5 = 20.0 25.0 2.3129773360E+01 1.7439951146E-01 b6 = 70.0 71.0 7.1994503004E+01 6.2622793913E-01 b7 = 178.0 180.0 1.7899805021E+02 1.2436988217E-01 b8 = 16.5 20.0 1.8389389025E+01 2.0134312832E-01 Residual Sum of Squares: 1.3158222432E+03 Residual Standard Deviation: 2.3317980180E+00 Degrees of Freedom: 242 Number of Observations: 250 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i250/; table data(i,*) y x i1 97.62227 1.000000 i2 97.80724 2.000000 i3 96.62247 3.000000 i4 92.59022 4.000000 i5 91.23869 5.000000 i6 95.32704 6.000000 i7 90.35040 7.000000 i8 89.46235 8.000000 i9 91.72520 9.000000 i10 89.86916 10.000000 i11 86.88076 11.00000 i12 85.94360 12.00000 i13 87.60686 13.00000 i14 86.25839 14.00000 i15 80.74976 15.00000 i16 83.03551 16.00000 i17 88.25837 17.00000 i18 82.01316 18.00000 i19 82.74098 19.00000 i20 83.30034 20.00000 i21 81.27850 21.00000 i22 81.85506 22.00000 i23 80.75195 23.00000 i24 80.09573 24.00000 i25 81.07633 25.00000 i26 78.81542 26.00000 i27 78.38596 27.00000 i28 79.93386 28.00000 i29 79.48474 29.00000 i30 79.95942 30.00000 i31 76.10691 31.00000 i32 78.39830 32.00000 i33 81.43060 33.00000 i34 82.48867 34.00000 i35 81.65462 35.00000 i36 80.84323 36.00000 i37 88.68663 37.00000 i38 84.74438 38.00000 i39 86.83934 39.00000 i40 85.97739 40.00000 i41 91.28509 41.00000 i42 97.22411 42.00000 i43 93.51733 43.00000 i44 94.10159 44.00000 i45 101.91760 45.00000 i46 98.43134 46.00000 i47 110.4214 47.00000 i48 107.6628 48.00000 i49 111.7288 49.00000 i50 116.5115 50.00000 i51 120.7609 51.00000 i52 123.9553 52.00000 i53 124.2437 53.00000 i54 130.7996 54.00000 i55 133.2960 55.00000 i56 130.7788 56.00000 i57 132.0565 57.00000 i58 138.6584 58.00000 i59 142.9252 59.00000 i60 142.7215 60.00000 i61 144.1249 61.00000 i62 147.4377 62.00000 i63 148.2647 63.00000 i64 152.0519 64.00000 i65 147.3863 65.00000 i66 149.2074 66.00000 i67 148.9537 67.00000 i68 144.5876 68.00000 i69 148.1226 69.00000 i70 148.0144 70.00000 i71 143.8893 71.00000 i72 140.9088 72.00000 i73 143.4434 73.00000 i74 139.3938 74.00000 i75 135.9878 75.00000 i76 136.3927 76.00000 i77 126.7262 77.00000 i78 124.4487 78.00000 i79 122.8647 79.00000 i80 113.8557 80.00000 i81 113.7037 81.00000 i82 106.8407 82.00000 i83 107.0034 83.00000 i84 102.46290 84.00000 i85 96.09296 85.00000 i86 94.57555 86.00000 i87 86.98824 87.00000 i88 84.90154 88.00000 i89 81.18023 89.00000 i90 76.40117 90.00000 i91 67.09200 91.00000 i92 72.67155 92.00000 i93 68.10848 93.00000 i94 67.99088 94.00000 i95 63.34094 95.00000 i96 60.55253 96.00000 i97 56.18687 97.00000 i98 53.64482 98.00000 i99 53.70307 99.00000 i100 48.07893 100.00000 i101 42.21258 101.00000 i102 45.65181 102.00000 i103 41.69728 103.00000 i104 41.24946 104.00000 i105 39.21349 105.00000 i106 37.71696 106.0000 i107 36.68395 107.0000 i108 37.30393 108.0000 i109 37.43277 109.0000 i110 37.45012 110.0000 i111 32.64648 111.0000 i112 31.84347 112.0000 i113 31.39951 113.0000 i114 26.68912 114.0000 i115 32.25323 115.0000 i116 27.61008 116.0000 i117 33.58649 117.0000 i118 28.10714 118.0000 i119 30.26428 119.0000 i120 28.01648 120.0000 i121 29.11021 121.0000 i122 23.02099 122.0000 i123 25.65091 123.0000 i124 28.50295 124.0000 i125 25.23701 125.0000 i126 26.13828 126.0000 i127 33.53260 127.0000 i128 29.25195 128.0000 i129 27.09847 129.0000 i130 26.52999 130.0000 i131 25.52401 131.0000 i132 26.69218 132.0000 i133 24.55269 133.0000 i134 27.71763 134.0000 i135 25.20297 135.0000 i136 25.61483 136.0000 i137 25.06893 137.0000 i138 27.63930 138.0000 i139 24.94851 139.0000 i140 25.86806 140.0000 i141 22.48183 141.0000 i142 26.90045 142.0000 i143 25.39919 143.0000 i144 17.90614 144.0000 i145 23.76039 145.0000 i146 25.89689 146.0000 i147 27.64231 147.0000 i148 22.86101 148.0000 i149 26.47003 149.0000 i150 23.72888 150.0000 i151 27.54334 151.0000 i152 30.52683 152.0000 i153 28.07261 153.0000 i154 34.92815 154.0000 i155 28.29194 155.0000 i156 34.19161 156.0000 i157 35.41207 157.0000 i158 37.09336 158.0000 i159 40.98330 159.0000 i160 39.53923 160.0000 i161 47.80123 161.0000 i162 47.46305 162.0000 i163 51.04166 163.0000 i164 54.58065 164.0000 i165 57.53001 165.0000 i166 61.42089 166.0000 i167 62.79032 167.0000 i168 68.51455 168.0000 i169 70.23053 169.0000 i170 74.42776 170.0000 i171 76.59911 171.0000 i172 81.62053 172.0000 i173 83.42208 173.0000 i174 79.17451 174.0000 i175 88.56985 175.0000 i176 85.66525 176.0000 i177 86.55502 177.0000 i178 90.65907 178.0000 i179 84.27290 179.0000 i180 85.72220 180.0000 i181 83.10702 181.0000 i182 82.16884 182.0000 i183 80.42568 183.0000 i184 78.15692 184.0000 i185 79.79691 185.0000 i186 77.84378 186.0000 i187 74.50327 187.0000 i188 71.57289 188.0000 i189 65.88031 189.0000 i190 65.01385 190.0000 i191 60.19582 191.0000 i192 59.66726 192.0000 i193 52.95478 193.0000 i194 53.87792 194.0000 i195 44.91274 195.0000 i196 41.09909 196.0000 i197 41.68018 197.0000 i198 34.53379 198.0000 i199 34.86419 199.0000 i200 33.14787 200.0000 i201 29.58864 201.0000 i202 27.29462 202.0000 i203 21.91439 203.0000 i204 19.08159 204.0000 i205 24.90290 205.0000 i206 19.82341 206.0000 i207 16.75551 207.0000 i208 18.24558 208.0000 i209 17.23549 209.0000 i210 16.34934 210.0000 i211 13.71285 211.0000 i212 14.75676 212.0000 i213 13.97169 213.0000 i214 12.42867 214.0000 i215 14.35519 215.0000 i216 7.703309 216.0000 i217 10.234410 217.0000 i218 11.78315 218.0000 i219 13.87768 219.0000 i220 4.535700 220.0000 i221 10.059280 221.0000 i222 8.424824 222.0000 i223 10.533120 223.0000 i224 9.602255 224.0000 i225 7.877514 225.0000 i226 6.258121 226.0000 i227 8.899865 227.0000 i228 7.877754 228.0000 i229 12.51191 229.0000 i230 10.66205 230.0000 i231 6.035400 231.0000 i232 6.790655 232.0000 i233 8.783535 233.0000 i234 4.600288 234.0000 i235 8.400915 235.0000 i236 7.216561 236.0000 i237 10.017410 237.0000 i238 7.331278 238.0000 i239 6.527863 239.0000 i240 2.842001 240.0000 i241 10.325070 241.0000 i242 4.790995 242.0000 i243 8.377101 243.0000 i244 6.264445 244.0000 i245 2.706213 245.0000 i246 8.362329 246.0000 i247 8.983658 247.0000 i248 3.362571 248.0000 i249 1.182746 249.0000 i250 4.875359 250.0000 ; * * 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' /9.8778210871E+01/ cb2 'certified value for b2' /1.0497276517E-02/ cb3 'certified value for b3' /1.0048990633E+02/ cb4 'certified value for b4' /6.7481111276E+01/ cb5 'certified value for b5' /2.3129773360E+01/ cb6 'certified value for b6' /7.1994503004E+01/ cb7 'certified value for b7' /1.7899805021E+02/ cb8 'certified value for b8' /1.8389389025E+01/ ce1 'certified std err for b1 ' /5.7527312730E-01/ ce2 'certified std err for b2 ' /1.1406289017E-04/ ce3 'certified std err for b3 ' /5.8831775752E-01/ ce4 'certified std err for b4 ' /1.0460593412E-01/ ce5 'certified std err for b5 ' /1.7439951146E-01/ ce6 'certified std err for b6 ' /6.2622793913E-01/ ce7 'certified std err for b7 ' /1.2436988217E-01/ ce8 'certified std err for b8 ' /2.0134312832E-01/ ; *----------------------------------------------------------------------------- * 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' b8 'coefficient to estimate' ; equations fit(i) 'the non-linear model' obj 'objective' ; obj.. sse =n= 0; fit(i).. y(i) =e= b1*exp(-b2*x(i)) +b3*exp(-sqr(x(i)-b4)/sqr(b5)) +b6*exp(-sqr(x(i)-b7)/sqr(b8)); option nlp=nls; model nlfit /obj,fit/; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l=9.7000000000E+01; b2.l=9.0000000000E-03; b3.l=1.0000000000E+02; b4.l=6.5000000000E+01; b5.l=2.0000000000E+01; b6.l=7.0000000000E+01; b7.l=1.7800000000E+02; b8.l=1.6500000000E+01; solve nlfit minimizing sse using nlp; 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)+abs(b8.l-cb8))>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)+abs(b8.m-ce8))>0.0001) "Accuracy problem"; *----------------------------------------------------------------------------- * second set of initial values *----------------------------------------------------------------------------- b1.l=9.4000000000E+01; b2.l=1.0500000000E-02; b3.l=9.9000000000E+01; b4.l=6.3000000000E+01; b5.l=2.5000000000E+01; b6.l=7.1000000000E+01; b7.l=1.8000000000E+02; b8.l=2.0000000000E+01; solve nlfit minimizing sse using nlp; 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)+abs(b8.l-cb8))>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)+abs(b8.m-ce8))>0.0001) "Accuracy problem"; *----------------------------------------------------------------------------- * third set of initial values *----------------------------------------------------------------------------- b1.l=9.8778210871E+01; b2.l=1.0497276517E-02; b3.l=1.0048990633E+02; b4.l=6.7481111276E+01; b5.l=2.3129773360E+01; b6.l=7.1994503004E+01; b7.l=1.7899805021E+02; b8.l=1.8389389025E+01; solve nlfit minimizing sse using nlp; 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)+abs(b8.l-cb8))>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)+abs(b8.m-ce8))>0.0001) "Accuracy problem";