'******************************************************* '* Program to demonstrate multi-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-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION * '* * '* How many data points ? 10 * '* * '* How many dimensions ? 2 * '* * '* What is the fit for dimension 1 ? 2 * '* What is the fit for dimension 2 ? 1 * '* * '* Input the data as prompted: * '* * '* Y( 1) = ? 7 * '* X( 1, 1) = ? 1 * '* X( 1, 2) = ? 6 * '* * '* Y( 2) = ? 7 * '* X( 2, 1) = ? 6 * '* X( 2, 2) = ? 1 * '* * '* Y( 3) = ? 6 * '* X( 3, 1) = ? 3 * '* X( 3, 2) = ? 3 * '* * '* Y( 4) = ? 8 * '* X( 4, 1) = ? 2 * '* X( 4, 2) = ? 6 * '* * '* Y( 5) = ? 9 * '* X( 5, 1) = ? 1 * '* X( 5, 2) = ? 8 * '* * '* Y( 6) = ? 9 * '* X( 6, 1) = ? 7 * '* X( 6, 2) = ? 2 * '* * '* Y( 7) = ? 6 * '* X( 7, 1) = ? 3 * '* X( 7, 2) = ? 3 * '* * '* Y( 8) = ? 7 * '* X( 8, 1) = ? 3 * '* X( 8, 2) = ? 4 * '* * '* Y( 9) = ? 7 * '* X( 9, 1) = ? 4 * '* X( 9, 2) = ? 3 * '* * '* Y( 10) = ? 2 * '* X( 10, 1) = ? 0 * '* X( 10, 2) = ? 2 * '* * '* The calculated coefficients are: * '* * '* 1 0 * '* 2 .999999 * '* 3 0 * '* 4 .999999 * '* 5 0 * '* 6 -1E-006 * '* * '* Standard deviation: 0 * '******************************************************* DEFINT I-N DEFDBL A-H, O-Z CLS PRINT " MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION" PRINT INPUT " How many data points ? ", m PRINT INPUT " How many dimensions ? ", l PRINT FOR i = 1 TO l PRINT " What is the fit for dimension "; i; " "; : INPUT m(i) NEXT i n = 1 FOR i = 1 TO l n = n * (m(i) + 1) NEXT i IF m < n THEN m = n DIM X(m, l), y(m), Z(m, n), d(n), A(m, m), B(m, 2 * m), C(m, m) PRINT PRINT " Input the data as prompted:" PRINT FOR i = 1 TO m PRINT " Y("; i; ") = "; : INPUT y(i) FOR j = 1 TO l PRINT " X("; i; ","; j; ") = "; : INPUT X(i, j) NEXT j PRINT NEXT i 'Call coefficients generation subroutine GOSUB 1000 'Call regression subroutine GOSUB 1200 PRINT " The calculated coefficients are:" PRINT FOR i = 1 TO n PRINT " "; i; " "; INT(1000000 * d(i)) / 1000000 NEXT i 'Get standard deviation n = n - 1 GOSUB 1000 PRINT PRINT " Standard deviation: "; INT(1000000 * d) / 1000000 PRINT END '************************************************************* '* Coefficient matrix generation subroutine * '* for multiple non-linear regression. * '* --------------------------------------------------------- * '* Also calculates the standard deviation d, even though * '* there is some redundant computing. * '* The maximum number of dimensions is 9. * '* The input data set consists of m data sets of the form: * '* Y(i),X(i,1),X(i,2) ... X(i,l) * '* The number of dimensions is l. * '* The order of the fit to each dimension is M(j). * '* The result is an (m1+1)(m2+1)...(ml+1)+1 column by m row * '* matrix, Z. This matrix is arranged as follows * '* (Ex.:l=2,M(1)=2,M(2)=2): * '* 1 X1 X1*X1 X2 X2*X1 X2*X1*X1 X2*X2 X2*X2*X1 X2*X2*X1*X1 * '* This matrix should be dimensioned in the calling program * '* as should also the X(i,j) matrix of data values. * '************************************************************* 'Calculate the total number of dimensions 1000 n = 1 FOR i = 1 TO l n = n * (m(i) + 1) NEXT d = 0 FOR i = 1 TO m 'Branch according to dimension l (return if l > 9) IF l > 0 THEN GOTO 10 l = 0: RETURN 10 IF l <= 9 THEN GOTO 15 l = 0: RETURN 15 j = 0 IF l = 1 THEN GOSUB 40 IF l = 2 THEN GOSUB 50 IF l = 3 THEN GOSUB 60 IF l = 4 THEN GOSUB 70 IF l = 5 THEN GOSUB 80 IF l = 6 THEN GOSUB 90 IF l = 7 THEN GOSUB 100 IF l = 8 THEN GOSUB 110 IF l = 9 THEN GOSUB 120 y = 0 FOR k = 1 TO n y = y + d(k) * Z(i, k) NEXT k d = d + (y(i) - y) * (y(i) - y) NEXT i 'Calculate standard deviation (if m > n) 30 IF m - n > 0 THEN GOTO 35 d = 0: RETURN 35 d = d / (m - n) d = SQR(d) RETURN 'Subroutines used by subroutine 1000 40 B = 1 41 C = B FOR i1 = 0 TO m(1) j = j + 1: Z(i, j) = B: B = B * X(i, 1) NEXT i1 B = C RETURN 50 B = 1 51 C = B FOR i2 = 0 TO m(2) GOSUB 41 B = B * X(i, 2) NEXT i2 B = C RETURN 60 B = 1 61 C = B FOR i3 = 0 TO m(3) GOSUB 51 B = B * X(i, 3) NEXT i3 B = C RETURN 70 B = 1 71 C = B FOR i4 = 0 TO m(4) GOSUB 61 B = B * X(i, 4) NEXT i4 B = C RETURN 80 B = 1 81 C = B FOR i5 = 0 TO m(5) GOSUB 71 B = B * X(i, 5) NEXT i5 B = C RETURN 90 B = 1 91 C = B FOR i6 = 0 TO m(6) GOSUB 81 B = B * X(i, 6) NEXT i6 B = C RETURN 100 B = 1 101 C = B FOR i7 = 0 TO m(7) GOSUB 91 B = B * X(i, 7) NEXT i7 B = C RETURN 110 B = 1 111 C = B FOR i8 = 0 TO m(8) GOSUB 101 B = B * X(i, 8) NEXT i8 B = C RETURN 120 B = 1 121 C = B FOR i9 = 0 TO m(9) GOSUB 111 B = B * X(i, 9) NEXT i9 B = C RETURN '********************************************************** '* Least squares fitting subroutine, general purpose * '* subroutine for multidimensional, nonlinear regression. * '* ------------------------------------------------------ * '* The equation fitted has the form: * '* Y = D(1)X1 + D(2)X2 + ... + D(n)Xn * '* The coefficients are returned by the program in D(i). * '* The X(i) can be simple powers of x, or functions. * '* Note that the X(i) are assumed to be independent. * '* The measured responses are Y(i), there are m of them. * '* Y is a m row column vector, Z(i,j) is a m by n matrix. * '* m must be >= n+2. The subroutine inputs are m, n, Y(i) * '* and Z(i,j) previously calculated. The subroutine calls * '* several other matrix routines during the calculation. * '********************************************************** 1200 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 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 MLTNLREG.BAS