' ******************************************************************
' Khi-2 test
' ******************************************************************

' ------------------------------------------------------------------
' 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 Khi2_Conform(      Obs()   AS INTEGER, _
                       Calc()  AS DOUBLE,  _
                       N_estim AS INTEGER, _
                 BYREF Khi2    AS DOUBLE,  _
                 BYREF DoF     AS INTEGER)

  DIM AS INTEGER Lb, Ub, I
  DIM AS DOUBLE  Diff

  Lb = LBOUND(Obs) : Ub = UBOUND(Obs)

  IF LBOUND(Calc) <> Lb OR UBOUND(Calc) <> Ub THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  Khi2 = 0
  ErrCode = FOk

  FOR I = Lb TO Ub
    Diff = Obs(I) - Calc(I)
    Khi2 = Khi2 + Diff * Diff / Calc(I)
  NEXT I

  DoF = Ub - Lb - N_estim
END SUB

SUB Khi2_Indep(      Obs() AS INTEGER, _
               BYREF Khi2  AS DOUBLE,  _
               BYREF DoF   AS INTEGER)

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2

  Lb1 = LBOUND(Obs, 1) : Ub1 = UBOUND(Obs, 1)
  Lb2 = LBOUND(Obs, 2) : Ub2 = UBOUND(Obs, 2)
 
  DIM AS INTEGER SumLin(Lb1 TO Ub1), SumCol(Lb2 TO Ub2), Sum, I, J

  DIM AS DOUBLE  Prob, Calc, Diff

  Sum = 0

  FOR I = Lb1 TO Ub1
    SumLin(I) = 0
    FOR J = Lb2 TO Ub2
      SumLin(I) = SumLin(I) + Obs(I, J)
    NEXT J
    Sum = Sum + SumLin(I)
  NEXT I

  FOR J = Lb2 TO Ub2
    SumCol(J) = 0
    FOR I = Lb1 TO Ub1
      SumCol(J) = SumCol(J) + Obs(I, J)
    NEXT I
  NEXT J

  Khi2 = 0
  FOR I = Lb1 TO Ub1
    Prob = SumLin(I) / Sum
    FOR J = Lb2 TO Ub2 
      Calc = SumCol(J) * Prob
      Diff = Obs(I, J) - Calc
      Khi2 = Khi2 + Diff * Diff / Calc
    NEXT J
  NEXT I

  DoF = (Ub2 - Lb2) * (Ub1 - Lb1)
END SUB
