'******************************************************* '* Program to demonstrate multidimensional operation * '* of the multi-nonlinear regression subroutine with * '* iterative error reduction * '* --------------------------------------------------- * '* Ref.: BASIC Scientific Subroutines Vol. II, By * '* F.R. Ruckdeschel, Byte/McGRAW-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 0 * '* * '* Standard deviation: 0 * '* * '* Number of iterations: 4 * '* * '******************************************************* DEFINT I-N DEFDBL A-H, O-Z CLS PRINT " MULTIDIMENSIONAL 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), d1(n), y1(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 iteration supervisor GOSUB 2000 PRINT PRINT " The calculated coefficients are:" PRINT FOR i = 1 TO n PRINT " "; i; " "; INT(1000000 * d(i)) / 1000000 NEXT i PRINT PRINT " Standard deviation: "; INT(1000000 * d) / 1000000 PRINT PRINT " Number of iterations: "; l1 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 '******************************************************************** '* Multi-dimensional polynomial regression iteration subroutine * '* ---------------------------------------------------------------- * '* This routine supervises the calling of several other subroutines * '* in order to iteratively fit least squares polynomials in more * '* one dimension. * '* The routine repeatedly calculates improved coefficients until * '* the standard deviation is no longer reduced. The inputs to the * '* subroutine are the number of dimensions l, the degree of fit for * '* each dimension m(i), and the input data, x(i) and y(i). * '* The coefficients are returned in d(i), with the standard devia- * '* tion in d. Also returned is the number of iterations tried, l1. * '* y1(i), d1(i) and d1 are used respectively to save the original * '* values of y(i) and the current values of d(i) and d. * '******************************************************************** 2000 l1 = 0 'Save the y(i) FOR i = 1 TO m y1(i) = y(i) NEXT 'Zero d1(i) FOR i = 1 TO n d1(i) = 0 NEXT 'Set the initial standard deviation high d1 = 10000000 'Call coefficients subroutine 2050 GOSUB 1000 'Call regression subroutine GOSUB 1200 'Get standard deviation GOSUB 1000 'If standard deviation is decreasing, continue IF d1 > d THEN GOTO 2100 'Terminate iteration FOR i = 1 TO n d(i) = d1(i) NEXT ' Restore y(i) FOR i = 1 TO m y(i) = y1(i) NEXT 'Get the final standard deviation GOSUB 1000 RETURN 'Save the standard deviation 2100 d1 = d: l1 = l1 + 1 'Augment coefficient matrix FOR i = 1 TO n d(i) = d1(i) + d(i) d1(i) = d(i) NEXT 'Restore y(i) FOR i = 1 TO m y(i) = y1(i) NEXT 'Reduce y(i) according to the d(i) GOSUB 2150 'We now have a set of error values GOTO 2050 'End subroutine 2000 'Subroutine 2150 2150 FOR i = 1 TO m j = 0 IF l = 1 THEN GOSUB 2160 IF l = 2 THEN GOSUB 2170 IF l = 3 THEN GOSUB 2180 IF l = 4 THEN GOSUB 2190 IF l = 5 THEN GOSUB 2200 IF l = 6 THEN GOSUB 2210 IF l = 7 THEN GOSUB 2220 IF l = 8 THEN GOSUB 2230 IF l = 9 THEN GOSUB 2240 'Array generated for row i y = 0 FOR k = 1 TO n y = y + d(k) * z(i, k) NEXT k y(i) = y(i) - y NEXT i RETURN 'End subroutine 2150 2160 b = 1 2161 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 2170 b = 1 2171 c = b FOR i2 = 0 TO m(2) GOSUB 2161: b = b * x(i, 2) NEXT i2 b = c RETURN 2180 b = 1 2181 c = b FOR i3 = 0 TO m(3) GOSUB 2171: b = b * x(i, 3) NEXT i3 b = c RETURN 2190 b = 1 2191 c = b FOR i4 = 0 TO m(4) GOSUB 2181: b = b * x(i, 4) NEXT i4 b = c RETURN 2200 b = 1 2201 c = b FOR i5 = 0 TO m(5) GOSUB 2191: b = b * x(i, 5) NEXT i5 b = c RETURN 2210 b = 1 2211 c = b FOR i6 = 0 TO m(6) GOSUB 2201: b = b * x(i, 6) NEXT i6 b = c RETURN 2220 b = 1 2221 c = b FOR i7 = 0 TO m(7) GOSUB 2211: b = b * x(i, 7) NEXT i7 b = c RETURN 2230 b = 1 2231 c = b FOR i8 = 0 TO m(7) GOSUB 2221: b = b * x(i, 7) NEXT i8 b = c RETURN 2240 b = 1 FOR i9 = 0 TO m(9) GOSUB 2231: b = b * x(i, 9) NEXT i9 RETURN 'End regiter.bas