' ******************************************************************
' Correlation and principal component analysis
' ******************************************************************
' Example taken from:
' P. DAGNELIE, Analyse statistique a plusieurs variables,
' Presses Agronomiques de Gembloux, Belgique, 1982
' ******************************************************************

#INCLUDE "fbmath.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

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

SUB PrintMatrix(Title AS STRING, A() AS DOUBLE)
' ------------------------------------------------------------------
' Print matrix A
' ------------------------------------------------------------------

  DIM AS INTEGER I, J

  PRINT : PRINT Title : PRINT

  FOR I = 1 TO UBOUND(A, 1)
    FOR J = 1 TO UBOUND(A, 2)
      PRINT USING "########.####"; A(I, J);
    NEXT J
    PRINT
  NEXT I
END SUB

SUB PrintVector(Title AS STRING, B() AS DOUBLE)
' ------------------------------------------------------------------
' Print vector B
' ------------------------------------------------------------------

  DIM AS INTEGER I

  PRINT : PRINT Title : PRINT

  FOR I = 1 TO UBOUND(B)
    PRINT USING "########.####"; B(I)
  NEXT I
END SUB

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

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

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

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

' Display results
PRINT
PrintVector "Mean vector", M()
PrintMatrix "Variance-covariance matrix", V()
PrintMatrix "Correlation matrix", R()

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

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

' 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 principal factors
PrinFac Z(), C(), F()

' Display results
PRINT
PrintVector "Eigenvalues of correlation matrix", Lambda()
PrintMatrix "Eigenvectors (columns) of correlation matrix", C()
PrintMatrix "Correlations between factors (columns) and variables (lines)", Rc()
PrintMatrix "Principal factors", F()

SLEEP
