$ontext Linear Least Squares Regression NIST test data Erwin kalvelagen, dec 2004 Reference: http://www.itl.nist.gov/div898/strd/lls/lls.shtml Wampler, R. H. (1970). A Report of the Accuracy of Some Widely-Used Least Squares Computer Programs. Journal of the American Statistical Association, 65, pp. 549-565. Model: Polynomial Class 6 Parameters (B0,B1,...,B5) y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5) Certified Regression Statistics Standard Deviation Parameter Estimate of Estimate B0 1.00000000000000 0.000000000000000 B1 0.100000000000000 0.000000000000000 B2 0.100000000000000E-01 0.000000000000000 B3 0.100000000000000E-02 0.000000000000000 B4 0.100000000000000E-03 0.000000000000000 B5 0.100000000000000E-04 0.000000000000000 Residual Standard Deviation 0.000000000000000 R-Squared 1.00000000000000 Certified Analysis of Variance Table Source of Degrees of Sums of Mean Variation Freedom Squares Squares F Statistic Regression 5 6602.91858365167 1320.58371673033 Infinity Residual 15 0.000000000000000 0.000000000000000 $offtext set i 'cases' /i1*i21/; table data(i,*) y x i1 1.00000 0 i2 1.11111 1 i3 1.24992 2 i4 1.42753 3 i5 1.65984 4 i6 1.96875 5 i7 2.38336 6 i8 2.94117 7 i9 3.68928 8 i10 4.68559 9 i11 6.00000 10 i12 7.71561 11 i13 9.92992 12 i14 12.75603 13 i15 16.32384 14 i16 20.78125 15 i17 26.29536 16 i18 33.05367 17 i19 41.26528 18 i20 51.16209 19 i21 63.00000 20 ; set j /j0*j5/; set j1(j); j1(j)$(ord(j)>1) = yes; parameter v(j); v(j) = ord(j)-1; parameter x(i,j); x(i,'j0') = 1; x(i,j1) = power(data(i,'x'),v(j1)); display x; variables b(j) 'coefficients to estimate' sse 'sum of squared errors' ; equation fit(i) 'equation to fit' sumsq ; sumsq.. sse =n= 0; fit(i).. data(i,'y') =e= sum(j, b(j)*x(i,j)); option lp = ls; model leastsq /fit,sumsq/; solve leastsq using lp minimizing sse; option decimals=8; display b.l; parameter Bcert(j) / j0 1.00000000000000 j1 0.100000000000000 j2 0.100000000000000E-01 j3 0.100000000000000E-02 j4 0.100000000000000E-03 j5 0.100000000000000E-04 /; scalar err "Sum of squared errors in estimates"; err = sum(j, sqr(bcert(j)-b.l(j))); display err; abort$(err>0.0001) "Solution not accurate";