' ******************************************************************
' Two-way analysis of variance (several observations per sample)
' ******************************************************************

#INCLUDE "fbmath.bi"

#DEFINE NA    2     ' Number of modalities of factor A
#DEFINE NB    2     ' Number of modalities of factor B
#DEFINE N    12     ' Number of observations for each sample
#DEFINE Alpha 0.05  ' Significance level

DIM AS INTEGER I, J, K  ' Loop variables

' The samples are stored in a 3D array Z, such that
' Z(I, J, K) contains the K-th observation for the
' I-th modality of factor A and the J-th modality of factor B

#IF __FB_VERSION__ < "0.17"
DIM AS DOUBLE Z(1 TO NA, 1 TO NB, 1 TO N) = _
{{4.9, 2.9}, {2.7, 3.9}, {4.6, 3.3}, {5.9, 4.8}, {4.1, 3.5}, {7.2, 6.1}, _
 {2.1, 2.2}, {1.1, 2.9}, {5.0, 3.5}, {2.4, 4.4}, {2.1, 3.0}, {3.9, 5.6}, _
 {4.5, 6.9}, {4.0, 5.4}, {1.9, 3.6}, {4.8, 3.3}, {7.5, 5.8}, {4.4, 6.0}, _
 {2.4, 3.6}, {4.8, 3.9}, {5.5, 5.0}, {6.8, 2.2}, {3.1, 5.0}, {4.1, 4.7}}
#ELSE
DIM AS DOUBLE Z(1 TO NA, 1 TO NB, 1 TO N) = _
{{{4.9, 2.9, 2.7, 3.9, 4.6, 3.3, 5.9, 4.8, 4.1, 3.5, 7.2, 6.1},  _
  {2.1, 2.2, 1.1, 2.9, 5.0, 3.5, 2.4, 4.4, 2.1, 3.0, 3.9, 5.6}}, _
 {{4.5, 6.9, 4.0, 5.4, 1.9, 3.6, 4.8, 3.3, 7.5, 5.8, 4.4, 6.0},  _
  {2.4, 3.6, 4.8, 3.9, 5.5, 5.0, 6.8, 2.2, 3.1, 5.0, 4.1, 4.7}}}
#ENDIF

SUB GetSample(Z() AS DOUBLE, I AS INTEGER, J AS INTEGER, X() AS DOUBLE)
' Get sample (I,J) from array Z into vector X

  DIM K AS INTEGER

  FOR K = 1 TO N
    X(K) = Z(I, J, K)
  NEXT K
END SUB

DIM AS DOUBLE  X(1 TO N)            ' Sample
DIM AS DOUBLE  M(1 TO NA, 1 TO NB)  ' Means
DIM AS DOUBLE  S(1 TO NA, 1 TO NB)  ' Standard dev.

DIM AS DOUBLE  V(1 TO 4)    ' Variances (A, B, interaction, residual)
DIM AS INTEGER DoF(1 TO 4)  ' Degrees of freedom (A, B, interaction, residual)
DIM AS DOUBLE  F(1 TO 3)    ' Variance ratios (A, B, interaction)
DIM AS DOUBLE  Fc(1 TO 3)   ' Critical values
DIM AS DOUBLE  P            ' Probability

' Compute means and SD's
FOR I = 1 TO NA
  FOR J = 1 TO NB
    GetSample Z(), I, J, X()
    M(I,J) = Mean(X())
    S(I,J) = StDev(X(), M(I,J))
  NEXT J
NEXT I

' Compare means
AnOVa2 N, M(), S(), V(), F(), DoF()

' Compute critical values
P = 1 - Alpha
FOR I = 1 TO 3
  Fc(I) = InvSnedecor(DoF(I), DoF(4), P)
NEXT I

' Print results
PRINT "Two-way ANOVA"
PRINT
PRINT USING "Source        Variance    D.o.F.     F      F(p = #.###)"; Alpha
PRINT "--------------------------------------------------------"
PRINT USING "Factor A     ####.####    ####  ####.####  ####.####"; V(1); DoF(1); F(1); Fc(1)
PRINT USING "Factor B     ####.####    ####  ####.####  ####.####"; V(2); DoF(2); F(2); Fc(2)
PRINT USING "Interaction  ####.####    ####  ####.####  ####.####"; V(3); DoF(3); F(3); Fc(3)
PRINT USING "Residual     ####.####    ####"; V(4); DoF(4)

SLEEP
