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