'****************************************************
'* 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