share/amatrix/wilcoxon.mac

   1 /* copyright 2007 by Robert Dodier
   2  * I release this file under the terms of the GNU General Public License.
   3  */
   4
   5 wilcoxon (X, C) := block ([X1, X0, X10, R, n1, n0, Rsum, W],
   6
   7     X1 : X [C = 1, 1],
   8     X0 : X [C = 0, 1],
   9     n1 : nrows (X1),
  10     n0 : nrows (X0),
  11     if n1 = 0 or n0 = 0
  12         then 1/2
  13         else
  14            (X10 : append (elements (X1), elements (X0)),
  15             R : ranks (X10),
  16             Rsum : sum (R[i], i, 1, n1),
  17             W : (Rsum - n1 * (n1 + 1) / 2) / (n1 * n0)));
  18
  19 ranks (L) := block ([n : length (L), x, dx, runs, %ranks, L1, L2],
  20     makelist ([L[i], i], i, 1, n),
  21     sort (%%, lambda ([x, y], orderlessp (first (x), first (y)))),
  22     [L1 : map (first, %%), L2 : map (second, %%)],
  23     makelist ([L2[i], i], i, 1, n),
  24     sort (%%, lambda ([x, y], orderlessp (first (x), first (y)))),
  25     %ranks : map (second, %%),
  26     runs : find_runs_nontrivial (L1),
  27     for e in runs do block ([avg_rank],
  28         avg_rank : (1/second(e)) * sum (%ranks [L2[i]], i, first(e), first(e) + second(e) - 1),
  29         for i : first(e) thru first(e) + second(e) - 1 do %ranks [L2[i]] : avg_rank),
  30     %ranks);
  31
  32 find_runs (x) :=
  33     if x = [] then []
  34     else block ([dx, ii, di],
  35         dx : makelist (x[i] - x[i - 1], i, 2, length (x)),
  36         ii : cons (1, 1 + ev (sublist_indices (dx, lambda ([x], x # 0 and x # 0.0 and x # 0b0)))),
  37         di : endcons (length (x) + 1 - last (ii), makelist (ii[i + 1] - ii[i], i, 1, length (ii) - 1)),
  38         makelist ([ii[i], di[i]], i, 1, length (di)));
  39
  40 find_runs_nontrivial (x) :=
  41    (find_runs (x),
  42     sublist (%%, lambda ([e], second (e) > 1)));
  43
  44 permute (L, p) :=
  45     map
  46        (second,
  47         sort
  48            (map (lambda ([x, i], [i, x]), L, p),
  49             lambda ([a, b], orderlessp (first(a),  first(b)))));
  50
  51 inverse_permutation (p) :=
  52     map
  53        (second,
  54         sort
  55            (makelist ([p[i], i], i, 1, length (p)),
  56             lambda ([a, b], orderlessp (first(a), first(b)))));