$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 monthly averaged atmospheric pressure differences between Easter Island and Darwin, Australia. This difference drives the trade winds in the southern hemisphere. Fourier analysis of the data reveals 3 significant cycles. The annual cycle is the strongest, but cycles with periods of approximately 44 and 26 months are also present. These cycles correspond to the El Nino and the Southern Oscillation. Arguments to the SIN and COS functions are in radians. Reference: Kahaner, D., C. Moler, and S. Nash, (1989). Numerical Methods and Software. Englewood Cliffs, NJ: Prentice Hall, pp. 441-445. Data: 1 Response (y = atmospheric pressure) 1 Predictor (x = time) 168 Observations Average Level of Difficulty Observed Data Model: Miscellaneous Class 9 Parameters (b1 to b9) y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 ) + b5*cos( 2*pi*x/b4 ) + b6*sin( 2*pi*x/b4 ) + b8*cos( 2*pi*x/b7 ) + b9*sin( 2*pi*x/b7 ) + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 11.0 10.0 1.0510749193E+01 1.7488832467E-01 b2 = 3.0 3.0 3.0762128085E+00 2.4310052139E-01 b3 = 0.5 0.5 5.3280138227E-01 2.4354686618E-01 b4 = 40.0 44.0 4.4311088700E+01 9.4408025976E-01 b5 = -0.7 -1.5 -1.6231428586E+00 2.8078369611E-01 b6 = -1.3 0.5 5.2554493756E-01 4.8073701119E-01 b7 = 25.0 26.0 2.6887614440E+01 4.1612939130E-01 b8 = -0.3 -0.1 2.1232288488E-01 5.1460022911E-01 b9 = 1.4 1.5 1.4966870418E+00 2.5434468893E-01 Residual Sum of Squares: 7.8853978668E+02 Residual Standard Deviation: 2.2269642403E+00 Degrees of Freedom: 159 Number of Observations: 168 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i168/; table data(i,*) y x i1 12.90000 1.000000 i2 11.30000 2.000000 i3 10.60000 3.000000 i4 11.20000 4.000000 i5 10.90000 5.000000 i6 7.500000 6.000000 i7 7.700000 7.000000 i8 11.70000 8.000000 i9 12.90000 9.000000 i10 14.30000 10.000000 i11 10.90000 11.00000 i12 13.70000 12.00000 i13 17.10000 13.00000 i14 14.00000 14.00000 i15 15.30000 15.00000 i16 8.500000 16.00000 i17 5.700000 17.00000 i18 5.500000 18.00000 i19 7.600000 19.00000 i20 8.600000 20.00000 i21 7.300000 21.00000 i22 7.600000 22.00000 i23 12.70000 23.00000 i24 11.00000 24.00000 i25 12.70000 25.00000 i26 12.90000 26.00000 i27 13.00000 27.00000 i28 10.90000 28.00000 i29 10.400000 29.00000 i30 10.200000 30.00000 i31 8.000000 31.00000 i32 10.90000 32.00000 i33 13.60000 33.00000 i34 10.500000 34.00000 i35 9.200000 35.00000 i36 12.40000 36.00000 i37 12.70000 37.00000 i38 13.30000 38.00000 i39 10.100000 39.00000 i40 7.800000 40.00000 i41 4.800000 41.00000 i42 3.000000 42.00000 i43 2.500000 43.00000 i44 6.300000 44.00000 i45 9.700000 45.00000 i46 11.60000 46.00000 i47 8.600000 47.00000 i48 12.40000 48.00000 i49 10.500000 49.00000 i50 13.30000 50.00000 i51 10.400000 51.00000 i52 8.100000 52.00000 i53 3.700000 53.00000 i54 10.70000 54.00000 i55 5.100000 55.00000 i56 10.400000 56.00000 i57 10.90000 57.00000 i58 11.70000 58.00000 i59 11.40000 59.00000 i60 13.70000 60.00000 i61 14.10000 61.00000 i62 14.00000 62.00000 i63 12.50000 63.00000 i64 6.300000 64.00000 i65 9.600000 65.00000 i66 11.70000 66.00000 i67 5.000000 67.00000 i68 10.80000 68.00000 i69 12.70000 69.00000 i70 10.80000 70.00000 i71 11.80000 71.00000 i72 12.60000 72.00000 i73 15.70000 73.00000 i74 12.60000 74.00000 i75 14.80000 75.00000 i76 7.800000 76.00000 i77 7.100000 77.00000 i78 11.20000 78.00000 i79 8.100000 79.00000 i80 6.400000 80.00000 i81 5.200000 81.00000 i82 12.00000 82.00000 i83 10.200000 83.00000 i84 12.70000 84.00000 i85 10.200000 85.00000 i86 14.70000 86.00000 i87 12.20000 87.00000 i88 7.100000 88.00000 i89 5.700000 89.00000 i90 6.700000 90.00000 i91 3.900000 91.00000 i92 8.500000 92.00000 i93 8.300000 93.00000 i94 10.80000 94.00000 i95 16.70000 95.00000 i96 12.60000 96.00000 i97 12.50000 97.00000 i98 12.50000 98.00000 i99 9.800000 99.00000 i100 7.200000 100.00000 i101 4.100000 101.00000 i102 10.60000 102.00000 i103 10.100000 103.00000 i104 10.100000 104.00000 i105 11.90000 105.00000 i106 13.60000 106.0000 i107 16.30000 107.0000 i108 17.60000 108.0000 i109 15.50000 109.0000 i110 16.00000 110.0000 i111 15.20000 111.0000 i112 11.20000 112.0000 i113 14.30000 113.0000 i114 14.50000 114.0000 i115 8.500000 115.0000 i116 12.00000 116.0000 i117 12.70000 117.0000 i118 11.30000 118.0000 i119 14.50000 119.0000 i120 15.10000 120.0000 i121 10.400000 121.0000 i122 11.50000 122.0000 i123 13.40000 123.0000 i124 7.500000 124.0000 i125 0.6000000 125.0000 i126 0.3000000 126.0000 i127 5.500000 127.0000 i128 5.000000 128.0000 i129 4.600000 129.0000 i130 8.200000 130.0000 i131 9.900000 131.0000 i132 9.200000 132.0000 i133 12.50000 133.0000 i134 10.90000 134.0000 i135 9.900000 135.0000 i136 8.900000 136.0000 i137 7.600000 137.0000 i138 9.500000 138.0000 i139 8.400000 139.0000 i140 10.70000 140.0000 i141 13.60000 141.0000 i142 13.70000 142.0000 i143 13.70000 143.0000 i144 16.50000 144.0000 i145 16.80000 145.0000 i146 17.10000 146.0000 i147 15.40000 147.0000 i148 9.500000 148.0000 i149 6.100000 149.0000 i150 10.100000 150.0000 i151 9.300000 151.0000 i152 5.300000 152.0000 i153 11.20000 153.0000 i154 16.60000 154.0000 i155 15.60000 155.0000 i156 12.00000 156.0000 i157 11.50000 157.0000 i158 8.600000 158.0000 i159 13.80000 159.0000 i160 8.700000 160.0000 i161 8.600000 161.0000 i162 8.600000 162.0000 i163 8.700000 163.0000 i164 12.80000 164.0000 i165 13.20000 165.0000 i166 14.00000 166.0000 i167 13.40000 167.0000 i168 14.80000 168.0000 ; * * extract data * parameter x(i),y(i); x(i) = data(i,'x'); y(i) = data(i,'y'); * * certified values * * * certified values * scalars cb1 'certified value for b1' / 1.0510749193E+01/ cb2 'certified value for b2' / 3.0762128085E+00/ cb3 'certified value for b3' / 5.3280138227E-01/ cb4 'certified value for b4' / 4.4311088700E+01/ cb5 'certified value for b5' /-1.6231428586E+00/ cb6 'certified value for b6' / 5.2554493756E-01/ cb7 'certified value for b7' / 2.6887614440E+01/ cb8 'certified value for b8' / 2.1232288488E-01/ cb9 'certified value for b8' / 1.4966870418E+00/ ce1 'certified std err for b1 ' /1.7488832467E-01/ ce2 'certified std err for b2 ' /2.4310052139E-01/ ce3 'certified std err for b3 ' /2.4354686618E-01/ ce4 'certified std err for b4 ' /9.4408025976E-01/ ce5 'certified std err for b5 ' /2.8078369611E-01/ ce6 'certified std err for b6 ' /4.8073701119E-01/ ce7 'certified std err for b7 ' /4.1612939130E-01/ ce8 'certified std err for b8 ' /5.1460022911E-01/ ce9 'certified std err for b9 ' /2.5434468893E-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' b9 'coefficient to estimate' ; equations fit(i) 'the non-linear model' obj 'objective' ; obj.. sse =n= 0; fit(i).. y(i) =e= b1 + b2*cos( 2*pi*x(i)/12 ) + b3*sin( 2*pi*x(i)/12 ) + b5*cos( 2*pi*x(i)/b4 ) + b6*sin( 2*pi*x(i)/b4 ) + b8*cos( 2*pi*x(i)/b7 ) + b9*sin( 2*pi*x(i)/b7 ); option nlp=nls; model nlfit /obj,fit/; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 11.0; b2.l = 3.0; b3.l = 0.5; b4.l = 40.0; b5.l = -0.7; b6.l = -1.3; b7.l = 25.0; b8.l = -0.3; b9.l = 1.4; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l,b6.l,b7.l,b8.l,b9.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)+abs(b8.l-cb8)+abs(b9.l-cb9))>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)+abs(b9.m-ce9))>0.0001) "Accuracy problem"; *----------------------------------------------------------------------------- * second set of initial values *----------------------------------------------------------------------------- b1.l = 10.0; b2.l = 3.0; b3.l = 0.5; b4.l = 44.0; b5.l = -1.5; b6.l = 0.5; b7.l = 26.0; b8.l = -0.1; b9.l = 1.5; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l,b6.l,b7.l,b8.l,b9.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)+abs(b8.l-cb8)+abs(b9.l-cb9))>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)+abs(b9.m-ce9))>0.0001) "Accuracy problem";