' ******************************************************************
' This program performs a least squares fit of a straight line:
'                        Y = B(0) + B(1) * X
' ******************************************************************

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

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

' Number of points
CONST N = 5

' Data
DIM AS DOUBLE X(1 TO N) = {10,     20,     30,    40,     50}
DIM AS DOUBLE Y(1 TO N) = { 0.1865, 0.3616, 0.537, 0.7359, 0.9238}

' Fitted parameters and variance-covariance matrix
DIM AS DOUBLE B(0 TO 1), V(0 TO 1, 0 TO 1)

' 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 variable
DIM AS INTEGER I

' ******************************************************************
' 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 "Linear regression: Y = B(0) + B(1) * X"
  PRINT Line1

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

  PRINT Line1

  FOR I = 0 TO 1
    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           = ###.########"; SGN(B(1)) * 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

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

  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 "Linear regression", "X", "Y"

  SetClipping

  COLOR 12

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

  LINE (Xpixel(Xmin), Ypixel(B(0) + B(1) * Xmin)) - _
       (Xpixel(Xmax), Ypixel(B(0) + B(1) * Xmax))

  CloseGraph
END SUB

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

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

' Compute predicted Y values
FOR I = 1 TO N
  Ycalc(I) = B(0) + B(1) * X(I)
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(N - 2, 1 - 0.5 * Alpha)
F = InvSnedecor(1, N - 2, 1 - Alpha)

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

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

SLEEP
