' ******************************************************************
' Two-way analysis of variance
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE FOk 0
' No error

#DEFINE MatErrDim -3
' Non-compatible dimensions

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ******************************************************************

SUB AnOVa2(N     AS INTEGER, _
           M()   AS DOUBLE,  _
           S()   AS DOUBLE,  _
           V()   AS DOUBLE,  _
           F()   AS DOUBLE,  _
           DoF() AS INTEGER)
' ------------------------------------------------------------------
' Input parameters : N     = number of observations for each sample
'                    M()   = matrix of means 
'                           (factor A as lines, factor B as columns)
'                    S()   = matrix of standard deviations
' Output parameters: V()   = variances 
'                           (factor A, factor B, interaction, residual)
'                    F()   = variance ratios 
'                           (factor A, factor B, interaction)
'                    DoF() = degrees of freedom
'                           (factor A, factor B, interaction, residual)                    
' ------------------------------------------------------------------

  DIM AS INTEGER Lb1, Lb2, Ub1, Ub2, N1, N2, I, J, P

  Lb1 = LBOUND(M, 1) : Lb2 = LBOUND(M, 2)
  Ub1 = UBOUND(M, 1) : Ub2 = UBOUND(M, 2)

  N1 = Ub1 - Lb1 + 1
  N2 = Ub2 - Lb2 + 1

  IF N1 < 2 OR N2 < 2 OR _ 
     Lb1 <> LBOUND(S, 1) OR Lb2 <> LBOUND(S, 2) OR _
     Ub1 <> UBOUND(S, 1) OR Ub2 <> UBOUND(S, 2) _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  DIM AS DOUBLE ML(Lb1 TO Ub1)  ' Line means
  DIM AS DOUBLE MC(Lb2 TO Ub2)  ' Column means
  DIM AS DOUBLE Xbar            ' Global mean
  DIM AS DOUBLE D               ' Difference of means
  DIM AS DOUBLE Sum             ' Intermediate sum
  DIM AS DOUBLE SS(0 TO 3)      ' Sum of squares (factorial, A, B, inter.)

  ErrCode = FOk

  ' Line means
  FOR I = 1 TO N1 
    Sum = 0
    FOR J = 1 TO N2 
      Sum = Sum + M(I, J)
    NEXT J
    ML(I) = Sum / N2
  NEXT I

  ' Column means
  FOR J = 1 TO N2 
    Sum = 0
    FOR I = 1 TO N1 
      Sum = Sum + M(I, J)
    NEXT I
    MC(J) = Sum / N1
  NEXT J

  ' Global mean  
  Sum = 0
  FOR I = 1 TO N1 
    Sum = Sum + ML(I)
  NEXT I
  Xbar = Sum / N1

  ' Residual variance
  IF N = 1 THEN
    V(4) = 0
  ELSE
    Sum = 0
    FOR I = 1 TO N1 
      FOR J = 1 TO N2 
        Sum = Sum + S(I, J) * S(I, J)
      NEXT J
    NEXT I
    P = N1 * N2
    DoF(4) = P * (N - 1)
    V(4) = Sum / P
  END IF

  ' Factorial sum of squares
  Sum = 0
  FOR I = 1 TO N1 
    FOR J = 1 TO N2 
      D = M(I, J) - Xbar  
      Sum = Sum + D * D
    NEXT J
  NEXT I
  SS(0) = N * Sum

  ' Factorial variance (factor A)
  Sum = 0
  FOR I = 1 TO N1 
    D = ML(I) - Xbar  
    Sum = Sum + D * D
  NEXT I
  SS(1) = N2 * N * Sum
  DoF(1) = N1 - 1
  V(1) = SS(1) / DoF(1)

  ' Factorial variance (factor B)
  Sum = 0
  FOR J = 1 TO N2 
    D = MC(J) - Xbar
    Sum = Sum + D * D
  NEXT J
  SS(2) = N1 * N * Sum
  DoF(2) = N2 - 1
  V(2) = SS(2) / DoF(2)

  ' Factorial variance (interaction)
  SS(3) = SS(0) - SS(1) - SS(2)
  DoF(3) = DoF(1) * DoF(2)
  V(3) = SS(3) / DoF(3)

  ' Variance ratios
  IF N = 1 THEN
    F(1) = V(1) / V(3)
    F(2) = V(2) / V(3)
  ELSE
    F(1) = V(1) / V(4)
    F(2) = V(2) / V(4)
    F(3) = V(3) / V(4)
  END IF
END SUB
