'**************************************************** '* Program to demonstrate least squares * '* polynomial fitting subroutine * '* ------------------------------------------------ * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 * '* [BIBLI 01]. * '* ------------------------------------------------ * '* SAMPLE RUN: * '* * '* LEAST SQUARES POLYNOMIAL FITTING * '* * '* What is the order of the fit : 4 * '* What is the error reduction factor: 0 * '* How many data pooints are there : 11 * '* * '* Input the data points as prompted: * '* * '* 1 X, Y = ? 1,1 * '* 2 X, Y = ? 2,2 * '* 3 X, Y = ? 3,3 * '* 4 X, Y = ? 4,4 * '* 5 X, Y = ? 5,5 * '* 6 X, Y = ? 6,6 * '* 7 X, Y = ? 7,7 * '* 8 X, Y = ? 8,8 * '* 9 X, Y = ? 9,9 * '* 10 X, Y = ? 0,0 * '* 11 X, Y = ? 1,1 * '* * '* Coefficients are: * '* 0 -1E-006 * '* 1 1 * '* 2 -1E-006 * '* 3 0 * '* 4 -1E-006 * '* * '* Standard deviation = 0 * '* * '**************************************************** DEFINT I-N DEFDBL A-H, O-Z CLS PRINT " LEAST SQUARES POLYNOMIAL FITTING" PRINT INPUT " What is the order of the fit : ", m INPUT " What is the error reduction factor: ", e INPUT " How many data pooints are there : ", n DIM x(n), y(n), v(n), a(n), b(n), c(n), d(n), c2(n), e(n), f(n) PRINT PRINT " Input the data points as prompted:" PRINT FOR i = 1 TO n PRINT " "; i; " X , Y = "; : INPUT x(i), y(i) NEXT i PRINT GOSUB 1000 PRINT "Coefficients are:" FOR i = 0 TO l PRINT " "; i; " "; INT(1000000 * c(i)) / 1000000 NEXT i PRINT PRINT " Standard deviation: "; INT(1000000 * d) / 1000000 PRINT END '***************************************************************** '* Least squares polynomial fitting subroutine * '* ------------------------------------------------------------- * '* This program least squares fits a polynomial to input data. * '* Forsythe orthogonal polynomials are used in the fitting. * '* The number of data points is n. * '* The data is input to the subroutine in x(i), y(i) pairs. * '* The coefficients are returned in c(i), * '* the smoothed data is returned in v(i), * '* the order of the fit is specified by m. * '* The standard deviation of the fit is returned in d. * '* There are two options available by use of the parameter e: * '* 1. if e=0, the fit is to order m, * '* 2. if e>0, the order of fit increases towards m, but will * '* stop if the relative standard deviation does not decrease * '* by more than e between successive fits. * '* The order of the fit then obtained is l. * '***************************************************************** 1000 n1 = m + 1 v1 = 10000000 'Initialize the arrays FOR i = 1 TO n1 a(i) = 0: b(i) = 0: f(i) = 0 NEXT FOR i = 1 TO n v(i) = 0: d(i) = 0 NEXT d1 = SQR(n) w = d1 FOR i = 1 TO n e(i) = 1 / w NEXT f1 = d1: a1 = 0 FOR i = 1 TO n a1 = a1 + x(i) * e(i) * e(i) NEXT c1 = 0 FOR i = 1 TO n c1 = c1 + y(i) * e(i) NEXT b(1) = 1 / f1: f(1) = b(1) * c1 FOR i = 1 TO n v(i) = v(i) + e(i) * c1 NEXT m = 1 10 'Save latest results FOR i = 1 TO l c2(i) = c(i) NEXT l2 = l: v2 = v: f2 = f1: a2 = a1: f1 = 0 FOR i = 1 TO n b1 = e(i) e(i) = (x(i) - a2) * e(i) - f2 * d(i) d(i) = b1 f1 = f1 + e(i) * e(i) NEXT f1 = SQR(f1) FOR i = 1 TO n e(i) = e(i) / f1 NEXT a1 = 0 FOR i = 1 TO n a1 = a1 + x(i) * e(i) * e(i) NEXT c1 = 0 FOR i = 1 TO n c1 = c1 + e(i) * y(i) NEXT m = m + 1: i = 0: 15 l = m - i: b2 = b(l): d1 = 0 IF l > 1 THEN d1 = b(l - 1) d1 = d1 - a2 * b(l) - f2 * a(l) b(l) = d1 / f1: a(l) = b2: i = i + 1 IF i <> m THEN GOTO 15 FOR i = 1 TO n v(i) = v(i) + e(i) * c1 NEXT FOR i = 1 TO n1 f(i) = f(i) + b(i) * c1 c(i) = f(i) NEXT v = 0 FOR i = 1 TO n v = v + (v(i) - y(i)) * (v(i) - y(i)) NEXT 'Note the division is by the number of degrees of freedom v = SQR(v / (n - l - 1)) l = m IF e = 0 THEN GOTO 20 'Test for minimal improvement IF ABS(v1 - v) / v < e THEN GOTO 50 'If error is larger, quit IF e * v > e * v1 THEN GOTO 50 v1 = v 20 IF m = n1 THEN GOTO 30 GOTO 10 30 'Shift the c(i) down, so c(0) is the constant term FOR i = 1 TO l c(i - 1) = c(i) NEXT c(l) = 0 'l is the order of the polynomial fitted l = l - 1: d = v RETURN 50 'Aborted sequence, recover last values l = l2: v = v2 FOR i = 1 TO l c(i) = c2(i) NEXT GOTO 30 RETURN 'End of file lsqrpoly.bas