lib/misc/wx-mp-sr.c

   1 #include <config.h>
   2 #include "wx-mp-sr.h"
   3
   4 /*********************************************************************
   5 *
   6 * Calculate the exact level of significance for a
   7 * Wilcoxon Matched-Pair Signed-Ranks Test using the sample's
   8 * Sum of Ranks W and the sample size (i.e., number of pairs) N.
   9 * This whole routine can be run as a stand-alone program.
  10 *
  11 * Use:
  12 * WX-MP-SR W N
  13 *
  14 * Copyright 1996, Rob van Son
  15 *
  16 * This program is free software; you can redistribute it and/or modify
  17 * it under the terms of the GNU General Public License as published by
  18 * the Free Software Foundation; either version 2 of the License, or
  19 * (at your option) any later version.
  20 *
  21 * This program is distributed in the hope that it will be useful,
  22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
  23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  24 * GNU General Public License for more details.
  25 *
  26 * You should have received a copy of the GNU General Public License along
  27 * with this program; if not, write to the Free Software Foundation, Inc.,
  28 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  29 *
  30 * -------------------------------------------------------
  31 *                 Rob van Son
  32 * Institute of Phonetic Sciences & IFOTT
  33 * University of Amsterdam, Spuistraat 210
  34 * NL-1012VT Amsterdam, The Netherlands
  35 * Tel.: (+31) 205252196 Fax.: (+31) 205252197
  36 * Email: r.j.j.h.vanson@uva.nl
  37 * WWW page: http://www.fon.hum.uva.nl/rob
  38 * -------------------------------------------------------
  39 *
  40 * This is the actual routine that calculates the exact (two-tailed)
  41 * level of significance for the Wilcoxon Matched-Pairs Signed-Ranks
  42 * test. The inputs are the Sum of Ranks of either the positive of
  43 * negative samples (W) and the sample size (N).
  44 * The Level of significance is calculated by checking for each
  45 * possible outcome (2**N possibilities) whether the sum of ranks
  46 * is larger than or equal to the observed Sum of Ranks (W).
  47 *
  48 * NOTE: The execution-time scales like ~ N*2**N, i.e., N*pow(2, N),
  49 * which is more than exponential. Adding a single pair to the sample
  50 * (i.e., increase N by 1) will more than double the time needed to
  51 * complete the calculations (apart from an additive constant).
  52 * The execution-time of this program can easily outrun your
  53 * patience.
  54 *
  55 ***********************************************************************/
  56
  57 double LevelOfSignificanceWXMPSR(double Winput, long int N)
  58 {
  59   unsigned long int W, MaximalW, NumberOfPossibilities, CountLarger;
  60   unsigned long int i, RankSum, j;
  61   double p;
  62
  63   /* Determine Wmax, i.e., work with the largest Rank Sum */
  64   MaximalW = N*(N+1)/2;
  65   if(Winput < MaximalW/2)Winput = MaximalW - Winput;
  66   W = Winput;    /* Convert to long int */
  67   if(W != Winput)++W;  /* Increase to next full integer */
  68
  69   /* The total number of possible outcomes is 2**N  */
  70   NumberOfPossibilities = 1 << N;
  71
  72   /* Initialize and loop. The loop-interior will be run 2**N times. */
  73   CountLarger = 0;
  74   /* Generate all distributions of sign over ranks as bit-patterns (i). */
  75   for(i=0; i < NumberOfPossibilities; ++i)
  76   {
  77     RankSum = 0;
  78     /*
  79        Shift "sign" bits out of i to determine the Sum of Ranks (j).
  80     */
  81     for(j=0; j < N; ++j)
  82     {
  83       if((i >> j) & 1)RankSum += j + 1;
  84     };
  85     /*
  86     * Count the number of "samples" that have a Sum of Ranks larger than
  87     * or equal to the one found (i.e., >= W).
  88     */
  89     if(RankSum >= W)++CountLarger;
  90   };
  91   /*****************************************************************
  92   * The level of significance is the number of outcomes with a
  93   * sum of ranks equal to or larger than the one found (W)
  94   * divided by the total number of possible outcomes.
  95   * The level is doubled to get the two-tailed result.
  96   ******************************************************************/
  97   p = 2*((double)CountLarger) / ((double)NumberOfPossibilities);
  98
  99   return p;
 100 }
 101