' ******************************************************************
' Non-parametric tests
' ******************************************************************

' ------------------------------------------------------------------
' 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 Ranks(      N()  AS INTEGER, _
                X()  AS DOUBLE,  _
                Sr() AS DOUBLE,  _
          BYREF Corr AS DOUBLE)
' ------------------------------------------------------------------
' Compute ranks for non-parametric tests
' ------------------------------------------------------------------ 
' Sr() = sum of ranks for each sample
' Corr = correction for ties = Sum k(k^2 - 1)  k = nb of ties
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, Lb1, I, J, I1, J1, NI, NE
  DIM AS DOUBLE  Y, Z, R

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

  IF LBOUND(X, 2) <> Lb OR UBOUND(X, 2) <> Ub OR _
     LBOUND(Sr)   <> Lb OR UBOUND(Sr)   <> Ub THEN
       ErrCode = MatErrDim
  END IF
   
  Lb1 = LBOUND(X, 1)

  Corr = 0
  ErrCode = FOk

  FOR I = Lb TO Ub
    Sr(I) = 0
  NEXT I

  FOR I = Lb TO Ub
    FOR J = Lb1 TO N(I) + Lb1 - 1
      Y = X(J, I)
      NE = 0                       ' Nb of values equal to Y
      NI = 0                       ' Nb of values less than Y
      FOR I1 = Lb TO Ub
        FOR J1 = Lb1 TO N(I1) + Lb1 - 1
          Z = X(J1, I1)
          IF Z < Y THEN
            NI = NI + 1
          ELSEIF Z = Y THEN
            NE = NE + 1
          END IF
        NEXT J1
      NEXT I1
      R = NI + (NE + 1) / 2        ' Mean rank of Y
      Sr(I) = Sr(I) + R            ' Sum of ranks for sample I
      IF NE > 1 THEN
        Corr = Corr + NE * (NE * NE - 1)
      END IF
    NEXT J
  NEXT I
END SUB

SUB Mann_Whitney(      X1() AS DOUBLE,  _
                       X2() AS DOUBLE,  _
                 BYREF U    AS DOUBLE,  _
                 BYREF Eps  AS DOUBLE)
' ------------------------------------------------------------------ 
' Mann-Whitney test
' ------------------------------------------------------------------ 

  DIM AS INTEGER Lb1, Lb2, Ub1, Ub2, Nmax, I, K
  DIM AS DOUBLE  Som, Prod, U1, U2, MU, VU, Corr

  Lb1 = LBOUND(X1) : Lb2 = LBOUND(X2)
  Ub1 = UBOUND(X1) : Ub2 = UBOUND(X2)

  DIM AS INTEGER N(1 TO 2)

  N(1) = Ub1 - Lb1 + 1
  N(2) = Ub2 - Lb2 + 1

  IF N(1) > N(2) THEN Nmax = N(1) ELSE Nmax = N(2)

  DIM AS DOUBLE X(1 TO Nmax, 1 TO 2)

  FOR I = Lb1 TO Ub1  ' Copy X1() into first column of X()
    X(I, 1) = X1(I)
  NEXT I

  FOR I = Lb2 TO Ub2  ' Copy X2() into second column of X()
    X(I, 2) = X2(I)
  NEXT I

  DIM AS DOUBLE Sr(1 TO 2)

  Ranks N(), X(), Sr(), Corr

  Som = N(1) + N(2)
  Prod = N(1) * N(2)
  U1 = Prod + N(1) * (N(1) + 1) / 2 - Sr(1)
  U2 = Prod + N(2) * (N(2) + 1) / 2 - Sr(2)
  IF U1 > U2 THEN U = U2 ELSE U = U1
  MU = Prod / 2
  VU = Prod * ((Som + 1) - Corr / Som / (Som - 1)) / 12
  Eps = (U - MU) / SQR(VU)
END SUB

SUB Wilcoxon(      X()   AS DOUBLE,  _
                   Y()   AS DOUBLE,  _
                   Ndiff AS INTEGER, _
             BYREF T     AS DOUBLE,  _
             BYREF Eps   AS DOUBLE)
' ------------------------------------------------------------------ 
' Wilcoxon test
' ------------------------------------------------------------------ 

  DIM AS INTEGER Lb, Ub, J, J1, J2, N
  DIM AS DOUBLE  Diff, MT, VT, Corr

  Lb = LBOUND(X)
  Ub = UBOUND(X)

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

  DIM AS INTEGER ND(1 TO 2)
  DIM AS DOUBLE  D(Lb TO Ub, 1 TO 2), Sr(1 TO 2)

  J1 = 0 
  J2 = 0
  FOR J = Lb TO Ub
    Diff = X(J) - Y(J)
    IF Diff < 0 THEN
      J1 = J1 + 1
      D(J1, 1) = ABS(Diff)  ' Negative difference
    ELSEIF Diff > 0 THEN
      J2 = J2 + 1
      D(J2, 2) = Diff       ' Positive difference
    END IF
  NEXT J

  ND(1) = J1       ' Nb of negative differences
  ND(2) = J2       ' Nb of positive differences
  Ndiff = J1 + J2  ' Nb of non-null differences

  Ranks ND(), D(), Sr(), Corr

  IF Sr(1) > Sr(2) THEN T = Sr(2) ELSE T = Sr(1)
  N = Ub - Lb + 1
  MT = N * (N + 1) / 4
  VT = MT * (2 * N + 1) / 6 - Corr / 48
  Eps = (T - MT) / SQR(VT)
END SUB

SUB Kruskal_Wallis(      N() AS INTEGER, _
                         X() AS DOUBLE,  _
                   BYREF H   AS DOUBLE,  _
                   BYREF DoF AS INTEGER)
' ------------------------------------------------------------------ 
' Kruskal-Wallis test
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, NT
  DIM AS DOUBLE  S, Corr

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

  DIM AS DOUBLE Sr(Lb TO Ub)

  Ranks N(), X(), Sr(), Corr

  IF ErrCode <> FOk THEN EXIT SUB

  S = 0 
  NT = 0

  FOR I = Lb TO Ub
    S = S + Sr(I) * Sr(I) / N(I)
    NT = NT + N(I)
  NEXT I

  H = 12 * S / NT / (NT + 1) - 3 * (NT + 1)
  H = H / (1 - Corr / NT / (NT * NT - 1))
  DoF = Ub - Lb
END SUB
