' ******************************************************************
' This program performs a polynomial least squares fit
' ******************************************************************

#INCLUDE "fbmath.bi"
#INCLUDE "fbgfx.bi"

' ******************************************************************
' General declarations
' ******************************************************************

' Number of observations
CONST N = 12

' Degree of polynomial
CONST Deg = 2

' Data
DIM AS DOUBLE X(1 TO N) = _
{0,      _
 0.0136, _
 0.023,  _
 0.0352, _
 0.048,  _
 0.075,  _
 0.1067, _
 0.1374, _
 0.1734, _
 0.2139, _
 0.2594, _
 0.3113}

DIM AS DOUBLE Y(1 TO N) = _
{3.97, _
 4.03, _
 4.1,  _
 4.2,  _
 4.28, _
 4.47, _
 4.66, _
 4.83, _
 4.99, _
 5.12, _
 5.25, _
 5.37}

' Fitted parameters and variance-covariance matrix
' B() is DIM SHARED to be used by function PltFunc
DIM SHARED AS DOUBLE B(0 TO Deg)
DIM AS DOUBLE V(0 TO Deg, 0 TO Deg)

' Computed Y values
DIM AS DOUBLE Ycalc(1 TO N)

' Statistical tests
DIM AS TRegTest Test

' Significance level
CONST Alpha = 0.05

' Student's t and Snedecor's F
DIM AS DOUBLE T, F

' Loop variables
DIM AS INTEGER I, J

' ******************************************************************
' Subprograms
' ******************************************************************

SUB WriteResults(X() AS DOUBLE, Y() AS DOUBLE, Ycalc() AS DOUBLE, _
                 B() AS DOUBLE, V() AS DOUBLE, Test AS TRegTest,  _
                 T AS DOUBLE, F AS DOUBLE)
' ------------------------------------------------------------------
' Writes results to screen
' ------------------------------------------------------------------

  DIM AS STRING Line1 = STRING$(74, "-")
  DIM AS STRING Line2 = STRING$(74, "=")

  DIM AS INTEGER I
  DIM AS DOUBLE  R(1 TO N), SB, Sr

  PRINT Line2
  PRINT "Polynomial regression: Y = B(0) + B(1) * X + B(2) * X^2 + ..."
  PRINT Line1

  PRINT "Parameter   Est.value         Std.dev.         ";
  PRINT USING "##"; 100 * (1 - Alpha);
  PRINT "% Confidence Interval"

  PRINT Line1

  FOR I = 0 TO Deg
    SB = SQR(V(I,I))
    PRINT USING "B(#)"; I;
    PRINT USING "########.########"; B(I); SB; B(I) - T * SB;
    PRINT ";";
    PRINT USING "########.########"; B(I) + T * SB
  NEXT I

  PRINT Line1

  Sr = SQR(Test.Vr)

  PRINT USING "Number of observations           : n           = ###"; N
  PRINT USING "Residual error                   : s           = ###.########"; Sr
  PRINT USING "Coefficient of correlation       : r           = ###.########"; SQR(Test.R2)
  PRINT USING "Coefficient of determination     : r^2         = ###.########"; Test.R2
  PRINT USING "Adjusted coef. of determination  : r^2(a)      = ###.########"; Test.R2a
  PRINT USING "Variance ratio (explained/resid) : F(###, ###) = #######.####"; Test.Nu1; Test.Nu2; Test.F
  PRINT USING "Critical variance ratio          : F(p = #.##) = #######.####"; 1 - Alpha; F

  PRINT Line1

  PRINT "  i        Y obs.       Y calc.      Residual      Std.res."

  PRINT Line1

  FOR I = 1 TO N
    R(I) = Y(I) - Ycalc(I)
    PRINT USING "###"; I;
    PRINT USING "#########.####"; Y(I); Ycalc(I); R(I); R(I) / Sr
  NEXT I

  PRINT Line2
END SUB

FUNCTION PltFunc(X AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Function to be plotted
' ------------------------------------------------------------------

  PltFunc = Poly(X, B(), Deg)
END FUNCTION

SUB PlotGraph(X() AS DOUBLE, Y() AS DOUBLE, B() AS DOUBLE)
' ------------------------------------------------------------------
' Plots observed points and fitted curve
' ------------------------------------------------------------------

  DIM AS DOUBLE Xmin, Xmax, Xstep  ' Ox scale
  DIM AS DOUBLE Ymin, Ymax, Ystep  ' Oy scale

  IF NOT InitScreen THEN
    PRINT "Unable to set graphic mode!"
    EXIT SUB
  END IF

  InitGraph

  AutoScale X(), LinScale, Xmin, Xmax, Xstep
  AutoScale Y(), LinScale, Ymin, Ymax, Ystep

  SetOxScale LinScale, Xmin, Xmax, Xstep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 4
  PlotOyAxis Xmin, 4

  PlotGrid

  WriteTitles "Polynomial regression", "X", "Y"

  SetClipping

  COLOR 12

  PlotCurve X(), Y(), 1, 3

  PlotFunc @PltFunc

  CloseGraph
END SUB

' ******************************************************************
' Main program
' ******************************************************************

' Perform regression
PolFit X(), Y(), B(), V()

IF MathErr <> MatOk THEN
  PRINT "Error in polynomial regression!"
  SLEEP
END IF

' Compute predicted Y values
FOR I = 1 TO N
  Ycalc(I) = Poly(X(I), B(), Deg)
NEXT I

' Update variance-covariance matrix and compute statistical tests
RegTest Y(), Ycalc(), V(), Test

' Compute Student's t and Snedecor's F
T = InvStudent(Test.Nu2, 1 - 0.5 * Alpha)
F = InvSnedecor(Test.Nu1, Test.Nu2, 1 - Alpha)

' Write results
WriteResults X(), Y(), Ycalc(), B(), V(), Test, T, F

' Plot Y values
PlotGraph X(), Y(), B()

SLEEP
