$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 slightly-blended 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 = 96.0 98.0 9.9018328406E+01 5.3748766879E-01 b2 = 0.009 0.0105 1.0994945399E-02 1.3335306766E-04 b3 = 103.0 103.0 1.0188022528E+02 5.9217315772E-01 b4 = 106.0 105.0 1.0703095519E+02 1.5006798316E-01 b5 = 18.0 20.0 2.3578584029E+01 2.2695595067E-01 b6 = 72.0 73.0 7.2045589471E+01 6.1721965884E-01 b7 = 151.0 150.0 1.5327010194E+02 1.9466674341E-01 b8 = 18.0 20.0 1.9525972636E+01 2.6416549393E-01 Residual Sum of Squares: 1.2475282092E+03 Residual Standard Deviation: 2.2704790782E+00 Degrees of Freedom: 242 Number of Observations: 250 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i250/; table data(i,*) y x i1 97.58776 1.000000 i2 97.76344 2.000000 i3 96.56705 3.000000 i4 92.52037 4.000000 i5 91.15097 5.000000 i6 95.21728 6.000000 i7 90.21355 7.000000 i8 89.29235 8.000000 i9 91.51479 9.000000 i10 89.60966 10.000000 i11 86.56187 11.00000 i12 85.55316 12.00000 i13 87.13054 13.00000 i14 85.67940 14.00000 i15 80.04851 15.00000 i16 82.18925 16.00000 i17 87.24081 17.00000 i18 80.79407 18.00000 i19 81.28570 19.00000 i20 81.56940 20.00000 i21 79.22715 21.00000 i22 79.43275 22.00000 i23 77.90195 23.00000 i24 76.75468 24.00000 i25 77.17377 25.00000 i26 74.27348 26.00000 i27 73.11900 27.00000 i28 73.84826 28.00000 i29 72.47870 29.00000 i30 71.92292 30.00000 i31 66.92176 31.00000 i32 67.93835 32.00000 i33 69.56207 33.00000 i34 69.07066 34.00000 i35 66.53983 35.00000 i36 63.87883 36.00000 i37 69.71537 37.00000 i38 63.60588 38.00000 i39 63.37154 39.00000 i40 60.01835 40.00000 i41 62.67481 41.00000 i42 65.80666 42.00000 i43 59.14304 43.00000 i44 56.62951 44.00000 i45 61.21785 45.00000 i46 54.38790 46.00000 i47 62.93443 47.00000 i48 56.65144 48.00000 i49 57.13362 49.00000 i50 58.29689 50.00000 i51 58.91744 51.00000 i52 58.50172 52.00000 i53 55.22885 53.00000 i54 58.30375 54.00000 i55 57.43237 55.00000 i56 51.69407 56.00000 i57 49.93132 57.00000 i58 53.70760 58.00000 i59 55.39712 59.00000 i60 52.89709 60.00000 i61 52.31649 61.00000 i62 53.98720 62.00000 i63 53.54158 63.00000 i64 56.45046 64.00000 i65 51.32276 65.00000 i66 53.11676 66.00000 i67 53.28631 67.00000 i68 49.80555 68.00000 i69 54.69564 69.00000 i70 56.41627 70.00000 i71 54.59362 71.00000 i72 54.38520 72.00000 i73 60.15354 73.00000 i74 59.78773 74.00000 i75 60.49995 75.00000 i76 65.43885 76.00000 i77 60.70001 77.00000 i78 63.71865 78.00000 i79 67.77139 79.00000 i80 64.70934 80.00000 i81 70.78193 81.00000 i82 70.38651 82.00000 i83 77.22359 83.00000 i84 79.52665 84.00000 i85 80.13077 85.00000 i86 85.67823 86.00000 i87 85.20647 87.00000 i88 90.24548 88.00000 i89 93.61953 89.00000 i90 95.86509 90.00000 i91 93.46992 91.00000 i92 105.8137 92.00000 i93 107.8269 93.00000 i94 114.0607 94.00000 i95 115.5019 95.00000 i96 118.5110 96.00000 i97 119.6177 97.00000 i98 122.1940 98.00000 i99 126.9903 99.00000 i100 125.7005 100.00000 i101 123.7447 101.00000 i102 130.6543 102.00000 i103 129.7168 103.00000 i104 131.8240 104.00000 i105 131.8759 105.00000 i106 131.9994 106.0000 i107 132.1221 107.0000 i108 133.4414 108.0000 i109 133.8252 109.0000 i110 133.6695 110.0000 i111 128.2851 111.0000 i112 126.5182 112.0000 i113 124.7550 113.0000 i114 118.4016 114.0000 i115 122.0334 115.0000 i116 115.2059 116.0000 i117 118.7856 117.0000 i118 110.7387 118.0000 i119 110.2003 119.0000 i120 105.17290 120.0000 i121 103.44720 121.0000 i122 94.54280 122.0000 i123 94.40526 123.0000 i124 94.57964 124.0000 i125 88.76605 125.0000 i126 87.28747 126.0000 i127 92.50443 127.0000 i128 86.27997 128.0000 i129 82.44307 129.0000 i130 80.47367 130.0000 i131 78.36608 131.0000 i132 78.74307 132.0000 i133 76.12786 133.0000 i134 79.13108 134.0000 i135 76.76062 135.0000 i136 77.60769 136.0000 i137 77.76633 137.0000 i138 81.28220 138.0000 i139 79.74307 139.0000 i140 81.97964 140.0000 i141 80.02952 141.0000 i142 85.95232 142.0000 i143 85.96838 143.0000 i144 79.94789 144.0000 i145 87.17023 145.0000 i146 90.50992 146.0000 i147 93.23373 147.0000 i148 89.14803 148.0000 i149 93.11492 149.0000 i150 90.34337 150.0000 i151 93.69421 151.0000 i152 95.74256 152.0000 i153 91.85105 153.0000 i154 96.74503 154.0000 i155 87.60996 155.0000 i156 90.47012 156.0000 i157 88.11690 157.0000 i158 85.70673 158.0000 i159 85.01361 159.0000 i160 78.53040 160.0000 i161 81.34148 161.0000 i162 75.19295 162.0000 i163 72.66115 163.0000 i164 69.85504 164.0000 i165 66.29476 165.0000 i166 63.58502 166.0000 i167 58.33847 167.0000 i168 57.50766 168.0000 i169 52.80498 169.0000 i170 50.79319 170.0000 i171 47.03490 171.0000 i172 46.47090 172.0000 i173 43.09016 173.0000 i174 34.11531 174.0000 i175 39.28235 175.0000 i176 32.68386 176.0000 i177 30.44056 177.0000 i178 31.98932 178.0000 i179 23.63330 179.0000 i180 23.69643 180.0000 i181 20.26812 181.0000 i182 19.07074 182.0000 i183 17.59544 183.0000 i184 16.08785 184.0000 i185 18.94267 185.0000 i186 18.61354 186.0000 i187 17.25800 187.0000 i188 16.62285 188.0000 i189 13.48367 189.0000 i190 15.37647 190.0000 i191 13.47208 191.0000 i192 15.96188 192.0000 i193 12.32547 193.0000 i194 16.33880 194.0000 i195 10.438330 195.0000 i196 9.628715 196.0000 i197 13.12268 197.0000 i198 8.772417 198.0000 i199 11.76143 199.0000 i200 12.55020 200.0000 i201 11.33108 201.0000 i202 11.20493 202.0000 i203 7.816916 203.0000 i204 6.800675 204.0000 i205 14.26581 205.0000 i206 10.66285 206.0000 i207 8.911574 207.0000 i208 11.56733 208.0000 i209 11.58207 209.0000 i210 11.59071 210.0000 i211 9.730134 211.0000 i212 11.44237 212.0000 i213 11.22912 213.0000 i214 10.172130 214.0000 i215 12.50905 215.0000 i216 6.201493 216.0000 i217 9.019605 217.0000 i218 10.80607 218.0000 i219 13.09625 219.0000 i220 3.914271 220.0000 i221 9.567886 221.0000 i222 8.038448 222.0000 i223 10.231040 223.0000 i224 9.367410 224.0000 i225 7.695971 225.0000 i226 6.118575 226.0000 i227 8.793207 227.0000 i228 7.796692 228.0000 i229 12.45065 229.0000 i230 10.61601 230.0000 i231 6.001003 231.0000 i232 6.765098 232.0000 i233 8.764653 233.0000 i234 4.586418 234.0000 i235 8.390783 235.0000 i236 7.209202 236.0000 i237 10.012090 237.0000 i238 7.327461 238.0000 i239 6.525136 239.0000 i240 2.840065 240.0000 i241 10.323710 241.0000 i242 4.790035 242.0000 i243 8.376431 243.0000 i244 6.263980 244.0000 i245 2.705892 245.0000 i246 8.362109 246.0000 i247 8.983507 247.0000 i248 3.362469 248.0000 i249 1.182678 249.0000 i250 4.875312 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.9018328406E+01/ cb2 'certified value for b2' /1.0994945399E-02/ cb3 'certified value for b3' /1.0188022528E+02/ cb4 'certified value for b4' /1.0703095519E+02/ cb5 'certified value for b5' /2.3578584029E+01/ cb6 'certified value for b6' /7.2045589471E+01/ cb7 'certified value for b7' /1.5327010194E+02/ cb8 'certified value for b8' /1.9525972636E+01/ ce1 'certified std err for b1 ' /5.3748766879E-01/ ce2 'certified std err for b2 ' /1.3335306766E-04/ ce3 'certified std err for b3 ' /5.9217315772E-01/ ce4 'certified std err for b4 ' /1.5006798316E-01/ ce5 'certified std err for b5 ' /2.2695595067E-01/ ce6 'certified std err for b6 ' /6.1721965884E-01/ ce7 'certified std err for b7 ' /1.9466674341E-01/ ce8 'certified std err for b8 ' /2.6416549393E-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.6000000000E+01; b2.l=9.0000000000E-03; b3.l=1.0300000000E+02; b4.l=1.0600000000E+02; b5.l=1.8000000000E+01; b6.l=7.2000000000E+01; b7.l=1.5100000000E+02; b8.l=1.8000000000E+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.8000000000E+01; b2.l=1.0500000000E-02; b3.l=1.0300000000E+02; b4.l=1.0500000000E+02; b5.l=2.0000000000E+01; b6.l=7.3000000000E+01; b7.l=1.5000000000E+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.9018328406E+01; b2.l=1.0994945399E-02; b3.l=1.0188022528E+02; b4.l=1.0703095519E+02; b5.l=2.3578584029E+01; b6.l=7.2045589471E+01; b7.l=1.5327010194E+02; b8.l=1.9525972636E+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";