$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 strongly-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 Average 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 = 94.9 96.0 9.8940368970E+01 5.3005192833E-01 b2 = 0.009 0.0096 1.0945879335E-02 1.2554058911E-04 b3 = 90.1 80.0 1.0069553078E+02 8.1256587317E-01 b4 = 113.0 110.0 1.1163619459E+02 3.5317859757E-01 b5 = 20.0 25.0 2.3300500029E+01 3.6584783023E-01 b6 = 73.8 74.0 7.3705031418E+01 1.2091239082E+00 b7 = 140.0 139.0 1.4776164251E+02 4.0488183351E-01 b8 = 20.0 25.0 1.9668221230E+01 3.7806634336E-01 Residual Sum of Squares: 1.2444846360E+03 Residual Standard Deviation: 2.2677077625E+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.60965 10.000000 i11 86.56187 11.00000 i12 85.55315 12.00000 i13 87.13053 13.00000 i14 85.67938 14.00000 i15 80.04849 15.00000 i16 82.18922 16.00000 i17 87.24078 17.00000 i18 80.79401 18.00000 i19 81.28564 19.00000 i20 81.56932 20.00000 i21 79.22703 21.00000 i22 79.43259 22.00000 i23 77.90174 23.00000 i24 76.75438 24.00000 i25 77.17338 25.00000 i26 74.27296 26.00000 i27 73.11830 27.00000 i28 73.84732 28.00000 i29 72.47746 29.00000 i30 71.92128 30.00000 i31 66.91962 31.00000 i32 67.93554 32.00000 i33 69.55841 33.00000 i34 69.06592 34.00000 i35 66.53371 35.00000 i36 63.87094 36.00000 i37 69.70526 37.00000 i38 63.59295 38.00000 i39 63.35509 39.00000 i40 59.99747 40.00000 i41 62.64843 41.00000 i42 65.77345 42.00000 i43 59.10141 43.00000 i44 56.57750 44.00000 i45 61.15313 45.00000 i46 54.30767 46.00000 i47 62.83535 47.00000 i48 56.52957 48.00000 i49 56.98427 49.00000 i50 58.11459 50.00000 i51 58.69576 51.00000 i52 58.23322 52.00000 i53 54.90490 53.00000 i54 57.91442 54.00000 i55 56.96629 55.00000 i56 51.13831 56.00000 i57 49.27123 57.00000 i58 52.92668 58.00000 i59 54.47693 59.00000 i60 51.81710 60.00000 i61 51.05401 61.00000 i62 52.51731 62.00000 i63 51.83710 63.00000 i64 54.48196 64.00000 i65 49.05859 65.00000 i66 50.52315 66.00000 i67 50.32755 67.00000 i68 46.44419 68.00000 i69 50.89281 69.00000 i70 52.13203 70.00000 i71 49.78741 71.00000 i72 49.01637 72.00000 i73 54.18198 73.00000 i74 53.17456 74.00000 i75 53.20827 75.00000 i76 57.43459 76.00000 i77 51.95282 77.00000 i78 54.20282 78.00000 i79 57.46687 79.00000 i80 53.60268 80.00000 i81 58.86728 81.00000 i82 57.66652 82.00000 i83 63.71034 83.00000 i84 65.24244 84.00000 i85 65.10878 85.00000 i86 69.96313 86.00000 i87 68.85475 87.00000 i88 73.32574 88.00000 i89 76.21241 89.00000 i90 78.06311 90.00000 i91 75.37701 91.00000 i92 87.54449 92.00000 i93 89.50588 93.00000 i94 95.82098 94.00000 i95 97.48390 95.00000 i96 100.86070 96.00000 i97 102.48510 97.00000 i98 105.7311 98.00000 i99 111.3489 99.00000 i100 111.0305 100.00000 i101 110.1920 101.00000 i102 118.3581 102.00000 i103 118.8086 103.00000 i104 122.4249 104.00000 i105 124.0953 105.00000 i106 125.9337 106.0000 i107 127.8533 107.0000 i108 131.0361 108.0000 i109 133.3343 109.0000 i110 135.1278 110.0000 i111 131.7113 111.0000 i112 131.9151 112.0000 i113 132.1107 113.0000 i114 127.6898 114.0000 i115 133.2148 115.0000 i116 128.2296 116.0000 i117 133.5902 117.0000 i118 127.2539 118.0000 i119 128.3482 119.0000 i120 124.8694 120.0000 i121 124.6031 121.0000 i122 117.0648 122.0000 i123 118.1966 123.0000 i124 119.5408 124.0000 i125 114.7946 125.0000 i126 114.2780 126.0000 i127 120.3484 127.0000 i128 114.8647 128.0000 i129 111.6514 129.0000 i130 110.1826 130.0000 i131 108.4461 131.0000 i132 109.0571 132.0000 i133 106.5308 133.0000 i134 109.4691 134.0000 i135 106.8709 135.0000 i136 107.3192 136.0000 i137 106.9000 137.0000 i138 109.6526 138.0000 i139 107.1602 139.0000 i140 108.2509 140.0000 i141 104.96310 141.0000 i142 109.3601 142.0000 i143 107.6696 143.0000 i144 99.77286 144.0000 i145 104.96440 145.0000 i146 106.1376 146.0000 i147 106.5816 147.0000 i148 100.12860 148.0000 i149 101.66910 149.0000 i150 96.44254 150.0000 i151 97.34169 151.0000 i152 96.97412 152.0000 i153 90.73460 153.0000 i154 93.37949 154.0000 i155 82.12331 155.0000 i156 83.01657 156.0000 i157 78.87360 157.0000 i158 74.86971 158.0000 i159 72.79341 159.0000 i160 65.14744 160.0000 i161 67.02127 161.0000 i162 60.16136 162.0000 i163 57.13996 163.0000 i164 54.05769 164.0000 i165 50.42265 165.0000 i166 47.82430 166.0000 i167 42.85748 167.0000 i168 42.45495 168.0000 i169 38.30808 169.0000 i170 36.95794 170.0000 i171 33.94543 171.0000 i172 34.19017 172.0000 i173 31.66097 173.0000 i174 23.56172 174.0000 i175 29.61143 175.0000 i176 23.88765 176.0000 i177 22.49812 177.0000 i178 24.86901 178.0000 i179 17.29481 179.0000 i180 18.09291 180.0000 i181 15.34813 181.0000 i182 14.77997 182.0000 i183 13.87832 183.0000 i184 12.88891 184.0000 i185 16.20763 185.0000 i186 16.29024 186.0000 i187 15.29712 187.0000 i188 14.97839 188.0000 i189 12.11330 189.0000 i190 14.24168 190.0000 i191 12.53824 191.0000 i192 15.19818 192.0000 i193 11.70478 193.0000 i194 15.83745 194.0000 i195 10.035850 195.0000 i196 9.307574 196.0000 i197 12.86800 197.0000 i198 8.571671 198.0000 i199 11.60415 199.0000 i200 12.42772 200.0000 i201 11.23627 201.0000 i202 11.13198 202.0000 i203 7.761117 203.0000 i204 6.758250 204.0000 i205 14.23375 205.0000 i206 10.63876 206.0000 i207 8.893581 207.0000 i208 11.55398 208.0000 i209 11.57221 209.0000 i210 11.58347 210.0000 i211 9.724857 211.0000 i212 11.43854 212.0000 i213 11.22636 213.0000 i214 10.170150 214.0000 i215 12.50765 215.0000 i216 6.200494 216.0000 i217 9.018902 217.0000 i218 10.80557 218.0000 i219 13.09591 219.0000 i220 3.914033 220.0000 i221 9.567723 221.0000 i222 8.038338 222.0000 i223 10.230960 223.0000 i224 9.367358 224.0000 i225 7.695937 225.0000 i226 6.118552 226.0000 i227 8.793192 227.0000 i228 7.796682 228.0000 i229 12.45064 229.0000 i230 10.61601 230.0000 i231 6.001000 231.0000 i232 6.765096 232.0000 i233 8.764652 233.0000 i234 4.586417 234.0000 i235 8.390782 235.0000 i236 7.209201 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.8940368970E+01/ cb2 'certified value for b2' /1.0945879335E-02/ cb3 'certified value for b3' /1.0069553078E+02/ cb4 'certified value for b4' /1.1163619459E+02/ cb5 'certified value for b5' /2.3300500029E+01/ cb6 'certified value for b6' /7.3705031418E+01/ cb7 'certified value for b7' /1.4776164251E+02/ cb8 'certified value for b8' /1.9668221230E+01/ ce1 'certified std err for b1 ' /5.3005192833E-01/ ce2 'certified std err for b2 ' /1.2554058911E-04/ ce3 'certified std err for b3 ' /8.1256587317E-01/ ce4 'certified std err for b4 ' /3.5317859757E-01/ ce5 'certified std err for b5 ' /3.6584783023E-01/ ce6 'certified std err for b6 ' /1.2091239082E+00/ ce7 'certified std err for b7 ' /4.0488183351E-01/ ce8 'certified std err for b8 ' /3.7806634336E-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= 94.9; b2.l= 0.009; b3.l= 90.1; b4.l= 113.0; b5.l= 20.0; b6.l= 73.8; b7.l= 140.0; b8.l= 20.0; 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= 96.0; b2.l= 0.0096; b3.l= 80.0; b4.l= 110.0; b5.l= 25.0; b6.l= 74.0; b7.l= 139.0; b8.l= 25.0; 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";