' ******************************************************************
' This program performs a multiple linear least squares fit
' ******************************************************************

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

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

' Number of observations
CONST N = 23

' Number of independent variables
CONST Nvar = 4

' Data
DIM AS DOUBLE X(1 TO N, 0 TO Nvar) = _
{{1, -0.27, 7.327, 0, 0}, _
 {1, -0.55, 7.4  , 0, 0}, _
 {1, -0.53, 7.74 , 0, 0}, _
 {1, -0.57, 7.95 , 0, 0}, _
 {1, -0.87, 7.9  , 0, 0}, _
 {1, -1.36, 7.931, 0, 0}, _
 {1, -0.39, 6.849, 0, 0}, _
 {1, -0.66, 7.508, 0, 0}, _
 {1, -0.33, 7.419, 0, 0}, _
 {1, -1.7 , 7.496, 0, 0}, _
 {1, -0.68, 7.027, 0, 0}, _
 {1, -0.79, 8.15 , 0, 0}, _
 {1, -0.82, 8.822, 0, 0}, _
 {1, -0.66, 8.334, 0, 0}, _
 {1,  0.02, 7.421, 0, 0}, _
 {1,  0.06, 7.862, 1, 0}, _
 {1, -0.3 , 8.483, 1, 0}, _
 {1,  0.07, 9.82 , 0, 0}, _
 {1,  0   , 7.641, 1, 0}, _
 {1, -0.8 , 7.601, 0, 1}, _
 {1, -1.05, 7.565, 0, 1}, _
 {1, -0.35, 7.993, 0, 0}, _
 {1, -0.11, 7.13 , 0, 0}}

DIM AS DOUBLE Y(1 TO N) = _
{3.21, _
 3.94, _
 3.66, _
 3.99, _
 4.06, _
 4.09, _
 3.36, _
 3.92, _
 3.58, _
 4.26, _
 3.06, _
 4.13, _
 4.27, _
 4.36, _
 3.72, _
 3.89, _
 4.39, _
 3.92, _
 3.89, _
 5.1 , _
 5.14, _
 3.68, _
 3.7 }

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

' 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 "Multilinear regression: Y = B(0) + B(1) * X1 + B(2) * X2 + ..."
  PRINT Line1

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

  PRINT Line1

  FOR I = 0 TO Nvar
    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

SUB PlotGraph(Y() AS DOUBLE, Ycalc() AS DOUBLE)
' ------------------------------------------------------------------
' Plots Y values (observed and estimated) and the
' theoretical line: Y calc. = Y obs.
' ------------------------------------------------------------------

  DIM AS DOUBLE Ymin, Ymax, Ystep  ' Common axis scale

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

  InitGraph

  AutoScale Y(), LinScale, Ymin, Ymax, Ystep

  SetOxScale LinScale, Ymin, Ymax, Ystep
  SetOyScale LinScale, Ymin, Ymax, Ystep

  PlotOxAxis Ymin, 4
  PlotOyAxis Ymin, 4

  PlotGrid

  WriteTitles "Multilinear regression", "Y obs.", "Y calc."

  SetClipping

  COLOR 12

  PlotCurve Y(), Ycalc(), 1, 3

  LINE (Xpixel(Ymin), Ypixel(Ymin)) - (Xpixel(Ymax), Ypixel(Ymax))

  CloseGraph
END SUB

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

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

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

' Compute predicted Y values
FOR I = 1 TO N
  Ycalc(I) = 0
  FOR J = 0 TO Nvar
    Ycalc(I) = Ycalc(I) + B(J) * X(I, J)
  NEXT J
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 Y(), Ycalc()

SLEEP
