'****************************************************
'*  Parabolic Least Squares Demonstration Program   *
'* ------------------------------------------------ *
'* Reference: BASIC Scientific Subroutines, Vol. II *
'*   By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981   *
'*   [BIBLI 01].                                    *
'* ------------------------------------------------ *
'* SAMPLE RUN:                                      *
'*                                                  *
'* This program calculates a parabolic              *
'* least squares fit to a given data set.           *
'*                                                  *
'* INSTRUCTIONS                                     *
'* ------------                                     *
'*                                                  *
'* The number of data coordinates provided must be  *
'* greater than three. Otherwise, a divide by zero  *
'* error may result.                                *
'*                                                  *
'* Number of data points : 4                        *
'*                                                  *
'* There are two input options:                     *
'* 1. input coordinate pairs                        *
'* 2. first input the independant variable values,  *
'*    then input dependant ones.                    *
'*    Your choice (1 or 2): 1                       *
'*                                                  *
'* ? 1 1                                            *
'* ? 2 2                                            *
'* ? 3 3                                            *
'* ? 5 5.01                                         *
'*                                                  *
'* Fitted equation is:                              *
'*   Y = -4.573E-003 + 1.002571 X                   *
'*                                                  *
'* Standard deviation of fit: 2.9E-003              *
'*                                                  *
'****************************************************
DEFINT I-N
DEFDBL A-H, O-Z
10 CLS
PRINT
PRINT "LEAST SQUARES CURVE FIT ROUTINE"
PRINT
PRINT "This program calculates a parabolic"
PRINT "least squares fit to a given data set"
PRINT
PRINT "INSTRUCTIONS"
PRINT "------------"
PRINT
PRINT " The number of data coordinates provided"
PRINT " must be greater than three. Otherwise, a"
PRINT " divide by zero error may result."
PRINT
INPUT " Number of data points : ", n
IF n < 4 THEN GOTO 10
DIM x(n + 1), y(n + 1)
20 CLS : PRINT
PRINT " There are two input options:"
PRINT " 1. input coordinate pairs"
PRINT " 2. first input the independant variable values, then input dependant ones."
INPUT " Your choice (1 or 2) : ", iz
PRINT
IF iz < 1 THEN GOTO 20
IF iz > 2 THEN GOTO 20
IF iz = 1 THEN GOTO 30
FOR m = 0 TO n - 1
  PRINT m + 1; : INPUT x(m)
NEXT m
PRINT
FOR m = 0 TO n - 1
  PRINT m + 1; : INPUT y(m)
NEXT m
GOTO 50
30 FOR m = 0 TO n - 1
  PRINT m + 1; : INPUT x(m), y(m)
NEXT m
REM Call linear least squares subroutine
50 GOSUB 1000
PRINT
PRINT " Fitted equation is:"
PRINT "   Y = "; INT(1000000 * a) / 1000000; " ";
IF b > 0 THEN PRINT "+";
PRINT INT(1000000 * b) / 1000000; "X";
IF c > 0 THEN PRINT " +";
PRINT INT(1000000 * c) / 1000000; "X^2";
PRINT:PRINT
PRINT " Standard deviation of fit: "; INT(10000 * d) / 10000
END

'************************************************
'*      Parabolic least squares subroutine      *
'* -------------------------------------------- *
'* The input data set is X(m), Y(m). The number *
'* of data points is n (n>3). The returned para-*
'* meters are a,b,c, coefficients of equation   *
'* Y = a + b X + c X^2, d, standard deviation   *
'* of fit.                                      *
'************************************************
1000 a0 = 1: a1 = 0: a2 = 0: a3 = 0: a4 = 0
b0 = 0: b1 = 0: b2 = 0
FOR m = 0 TO n - 1
  a1 = a1 + x(m)
  a2 = a2 + x(m) * x(m)
  a3 = a3 + x(m) * x(m) * x(m)
  a4 = a4 + x(m) * x(m) * x(m) * x(m)
  b0 = b0 + y(m)
  b1 = b1 + y(m) * x(m)
  b2 = b2 + y(m) * x(m) * x(m)
NEXT m
a1 = a1 / n: a2 = a2 / n: a3 = a3 / n: a4 = a4 / n
b0 = b0 / n: b1 = b1 / n: b2 = b2 / n
d = a0 * (a2 * a4 - a3 * a3) - a1 * (a1 * a4 - a2 * a3) + a2 * (a1 * a3 - a2 * a2)
a = b0 * (a2 * a4 - a3 * a3) + b1 * (a2 * a3 - a1 * a4) + b2 * (a1 * a3 - a2 * a2)
a = a / d
b = b0 * (a2 * a3 - a1 * a4) + b1 * (a0 * a4 - a2 * a2) + b2 * (a1 * a2 - a0 * a3)
b = b / d
c = b0 * (a1 * a3 - a2 * a2) + b1 * (a2 * a1 - a0 * a3) + b2 * (a0 * a2 - a1 * a1)
c = c / d
'Evaluation of standard deviation d
d = 0
FOR m = 0 TO n - 1
  d1 = y(m) - a - b * x(m) - c * x(m) * x(m)
  d = d + d1 * d1
NEXT m
d = SQR(d / (n - 3))
RETURN

'End lstsqr2.bas