' ******************************************************************
' One-way analysis of variance
' ******************************************************************

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

#DEFINE FOk 0
' No error

#DEFINE FSing -2
' Singularity

#DEFINE MatErrDim -3
' Non-compatible dimensions

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

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

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

SUB AnOVa1(      N()   AS INTEGER, _
                 M()   AS DOUBLE,  _
                 S()   AS DOUBLE,  _
           BYREF V_f   AS DOUBLE,  _ 
           BYREF V_r   AS DOUBLE,  _ 
           BYREF F     AS DOUBLE,  _
           BYREF DoF_f AS INTEGER, _
           BYREF DoF_r AS INTEGER)
' ------------------------------------------------------------------
' Input parameters : N()          = samples sizes
'                    M()          = samples means
'                    S()          = samples SD's (computed with StDev) 
' Output parameters: V_f, V_r     = variances (factorial, residual)
'                    F            = ratio Vf / Vr 
'                    DoF_f, DoF_r = degrees of freedom
' ------------------------------------------------------------------    

  DIM AS INTEGER Lb, Ub, I, Ns, Nt

  Lb = LBOUND(N)
  Ub = UBOUND(N)

  Ns = Ub - Lb + 1  ' Number of samples

  IF Ns < 2 OR LBOUND(M) <> Lb OR UBOUND(M) <> Ub _
    OR LBOUND(S) <> Lb OR UBOUND(S) <> Ub _
  THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  Nt = 0
  FOR I = Lb TO Ub
    Nt = Nt + N(I)
  NEXT I

  IF Nt <= Ns THEN
    ErrCode = FSing
    EXIT SUB
  END IF

  DIM AS DOUBLE Xbar      ' Global mean
  DIM AS DOUBLE SSf, SSr  ' Sum of squares
  DIM AS DOUBLE D         ' Difference of means

  ErrCode = FOk

  Xbar = 0
  FOR I = Lb TO Ub
    Xbar = Xbar + N(I) * M(I)
  NEXT I
  Xbar = Xbar / Nt

  SSf = 0
  SSr = 0
  FOR I = Lb TO Ub
    D = M(I) - Xbar
    SSf = SSf + N(I) * D * D
    SSr = SSr + (N(I) - 1) * S(I) * S(I)
  NEXT I

  DoF_f = Ns - 1
  DoF_r = Nt - Ns
  V_f = SSf / DoF_f
  V_r = SSr / DoF_r
  F = V_f / V_r
END SUB
