'**************************************************** '* Program to demonstrate one dimensional operation * '* of the multi-nonlinear regression subroutine * '* ------------------------------------------------ * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 * '* [BIBLI 01]. * '* ------------------------------------------------ * '* SAMPLE RUN: * '* * '* MULTI-NON LINEAR REGRESSION * '* * '* Number of data points: 11 * '* Degree of the polynomial to be fitted: 4 * '* * '* Input the data in (x,y) pairs 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 * '* * '* The calculated coefficients are: * '* * '* 1 0 * '* 2 .999999 * '* 3 0 * '* 4 -1E-006 * '* 5 0 * '* * '* Standard deviation: 0 * '* * '**************************************************** DEFINT I-N DEFDBL A-H, O-Z CLS PRINT " MULTI-NONLINEAR REGRESSION" PRINT INPUT " Number of data points: ", m INPUT " Degree of the polynomial to be fitted: ", n DIM x(m), y(m), Z(m, n + 1), d(n + 1), A(m, m), B(m, 2 * m), C(m, m) PRINT PRINT " Input the data in (x,y) pairs as prompted:" PRINT FOR i = 1 TO m PRINT i; " X , Y = "; : INPUT x(i), y(i) NEXT PRINT 'Call coefficients generation subroutine GOSUB 1000 'Call regression subroutine n = n + 1 GOSUB 1200 PRINT " The calculated coefficients are:" FOR i = 1 TO n PRINT i; " "; INT(1000000 * d(i)) / 1000000 NEXT n = n - 1 'Call standard deviation subroutine GOSUB 1500 PRINT PRINT " Standard deviation: "; INT(1000000 * d) / 1000000 PRINT END 1000 'Coefficient matrix generation FOR i = 1 TO m B = 1 FOR j = 1 TO n + 1 Z(i, j) = B B = B * x(i) NEXT j NEXT i RETURN 1200 'Least squares multidimensional fitting subroutine m4 = m n4 = n FOR i = 1 TO m FOR j = 1 TO n A(i, j) = Z(i, j) NEXT j NEXT i GOSUB 5100 'b=Transpose(a) n1 = m: n2 = n: GOSUB 5400 'move A to C n1 = n: n2 = m: GOSUB 5200 'move B to A n1 = m: n2 = n: GOSUB 5300 'move C to B m1 = n: n1 = m: n2 = n: GOSUB 5000 'multiply A and B n1 = n: GOSUB 5500 'move C to A GOSUB 6000 'b=Inverse(a) m = m4 'restore m GOSUB 5200 'move B to A FOR i = 1 TO m FOR j = 1 TO n B(j, i) = Z(i, j) NEXT j NEXT i m2 = n: n2 = m: GOSUB 5000 'multiply A and B n1 = n: n2 = m: GOSUB 5500 'move C to A FOR i = 1 TO m B(i, 1) = y(i) NEXT m1 = n: n2 = 1: n1 = m: GOSUB 5000 'multiply A and B 'Product C is N by 1 - Regression coefficients are in C(I,1) FOR i = 1 TO n d(i) = C(i, 1) NEXT RETURN 1500 'Standard deviation calculation for a polynomial fit d = 0 FOR i = 1 TO m y = 0: B = 1 FOR j = 1 TO n + 1 y = y + d(j) * B B = B * x(i) NEXT j d = d + (y - y(i)) * (y - y(i)) NEXT i IF m - n - 1 > 0 THEN GOTO 1510 d = 0: RETURN 1510 d = d / (m - n - 1) d = SQR(d) RETURN 5000 'Matrix multiplication FOR i = 1 TO m1 FOR j = 1 TO n2 C(i, j) = 0 FOR k = 1 TO n1 C(i, j) = C(i, j) + A(i, k) * B(k, j) NEXT k NEXT j NEXT i RETURN 5100 'Matrix transpose FOR i = 1 TO n FOR j = 1 TO m B(i, j) = A(j, i) NEXT j NEXT i RETURN 5200 'Matrix save (B in A) IF n1 * n2 = 0 THEN RETURN FOR i1 = 1 TO n1 FOR i2 = 1 TO n2 A(i1, i2) = B(i1, i2) NEXT i2 NEXT i1 RETURN 5300 'Matrix save (C in B) IF n1 * n2 = 0 THEN RETURN FOR i1 = 1 TO n1 FOR i2 = 1 TO n2 B(i1, i2) = C(i1, i2) NEXT i2 NEXT i1 RETURN 5400 'Matrix save (A in C) IF n1 * n2 = 0 THEN RETURN FOR i1 = 1 TO n1 FOR i2 = 1 TO n2 C(i1, i2) = A(i1, i2) NEXT i2 NEXT i1 RETURN 5500 'Matrix save (C in A) IF n1 * n2 = 0 THEN RETURN FOR i1 = 1 TO n1 FOR i2 = 1 TO n2 A(i1, i2) = C(i1, i2) NEXT i2 NEXT i1 RETURN 6000 'Matrix inversion FOR i = 1 TO n FOR j = 1 TO n B(i, j + n) = 0 B(i, j) = A(i, j) NEXT j B(i, i + n) = 1 NEXT i FOR k = 1 TO n IF k = n THEN GOTO 6010 m = k FOR i = k + 1 TO n IF ABS(B(i, k)) > ABS(B(m, k)) THEN m = i NEXT i IF m = k THEN GOTO 6010 FOR j = k TO 2 * n B = B(k, j) B(k, j) = B(m, j) B(m, j) = B NEXT j 6010 FOR j = k + 1 TO 2 * n B(k, j) = B(k, j) / B(k, k) NEXT j IF k = 1 THEN GOTO 6020 FOR i = 1 TO k - 1 FOR j = k + 1 TO 2 * n B(i, j) = B(i, j) - B(i, k) * B(k, j) NEXT j NEXT i IF k = n THEN GOTO 6030 6020 FOR i = k + 1 TO n FOR j = k + 1 TO 2 * n B(i, j) = B(i, j) - B(i, k) * B(k, j) NEXT j NEXT i NEXT k 6030 FOR i = 1 TO n FOR j = 1 TO n B(i, j) = B(i, j + n) NEXT j NEXT i RETURN 'End Lstsqr.bas