' ******************************************************************
' Principal component analysis: scores plot
' ******************************************************************
' Example taken from:
' P. DAGNELIE, Analyse statistique a plusieurs variables,
' Presses Agronomiques de Gembloux, Belgique, 1982
' ******************************************************************

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

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

' Number of observations
CONST N = 11

' Number of independent variables
CONST Nvar = 4

' Data
DIM AS DOUBLE X(1 TO N, 1 TO Nvar) = _
{{ 87.9, 19.6,   1  , 1661}, _
 { 89.9, 15.2,  90.1,  968}, _
 {153  , 19.7,  56.6, 1353}, _
 {132.1, 17  ,  91  , 1293}, _
 { 88.8, 18.3,  93.7, 1153}, _
 {220.9, 17.8, 106.9, 1286}, _
 {117.7, 17.8,  65.5, 1104}, _
 {109  , 18.3,  41.8, 1574}, _
 {156.1, 17.8,  57.4, 1222}, _
 {181.5, 16.8, 140.6,  902}, _
 {181.4, 17  ,  74.3, 1150}}

DIM AS DOUBLE M(1 TO Nvar)              ' Mean vector
DIM AS DOUBLE V(1 TO Nvar, 1 TO Nvar)   ' Variance-covariance matrix
DIM AS DOUBLE R(1 TO Nvar, 1 TO Nvar)   ' Correlation matrix

DIM AS DOUBLE S(1 TO Nvar)              ' Standard deviations
DIM AS DOUBLE Lambda(1 TO Nvar)         ' Eigenvalues of correlation matrix

DIM AS DOUBLE C(1 TO Nvar, 1 TO Nvar)   ' Eigenvectors of correlation matrix
DIM AS DOUBLE Rc(1 TO Nvar, 1 TO Nvar)  ' Correlation factors/variables

DIM AS DOUBLE Z(1 TO N, 1 TO Nvar)      ' Scaled variables
DIM AS DOUBLE F(1 TO N, 1 TO Nvar)      ' Principal factors

' ******************************************************************
' Subprogram
' ******************************************************************

SUB PlotScores(F() AS DOUBLE)
' ------------------------------------------------------------------
' Plots the coordinates ('scores') of the original observations
' in the plane defined by the first two principal factors
' ------------------------------------------------------------------

  DIM AS DOUBLE  F1(1 TO N), F2(1 TO N)  ' Scores

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

  DIM AS INTEGER X0, Y0, X1, Y1, X2, Y2, I

  ' Extract scores from F matrix
  FOR I = 1 TO N
    F1(I) = F(I, 1)
    F2(I) = F(I, 2)
  NEXT I

  ' Determine scale
  AutoScale F1(), LinScale, Xmin, Xmax, Xstep
  AutoScale F2(), LinScale, Ymin, Ymax, Ystep

  ' Make sure the scales are centered on 0
  Xmax = Max(Xmax, ABS(Xmin)) : Xmin = -Xmax
  Ymax = Max(Ymax, ABS(Ymin)) : Ymin = -Ymax

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

  InitGraph

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

  PlotOxAxis Ymin, 2
  PlotOyAxis Xmin, 2

  PlotGrid BothGrid

  WriteTitles "PCA scores", "F1", "F2"

  X0 = Xpixel(0)
  Y0 = Ypixel(0)

  X1 = Xpixel(Xmin)
  Y1 = Ypixel(Ymin)

  X2 = Xpixel(Xmax)
  Y2 = Ypixel(Ymax)

  LINE (X1, Y0)-(X2, Y0), 1
  LINE (X0, Y1)-(X0, Y2), 1

  COLOR 12

  PlotCurve F1(), F2(), 1, 3

  CloseGraph
END SUB

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

' Compute mean vector
VecMean X(), M()

' Compute variance-covariance matrix
MatVarCov X(), M(), V()

' Compute correlation matrix
MatCorrel V(), R()

' Perform principal component analysis
' The original matrix R is destroyed
PCA R(), Lambda(), C(), Rc()

IF MathErr() = MatNonConv THEN
  PRINT "Non-convergence of eigenvalue computation"
  STOP
END IF

' Compute standard deviations
VecSD X(), M(), S()

' Scale variables
ScaleVar X(), M(), S(), Z()

' Compute principal factors
PrinFac Z(), C(), F()

' Plot graphic
PlotScores F()

