/* Wilcoxon rank sum test (two-sample nonparametric) */
/* Normal approximation of test statistic */
/* Hodges-Lehmann estimator of difference between sample means (medians) */
/* HL-like nonparametric confidence interval for same */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "statlib.h"

#define PROG_DATE   "January 31, 2007"
#define OUTPUT_NAME "twosamps_ranks.txt"

#define MAXMN 500        /* Enough for mm, nn, mm+nn, and mm*nn */

/* First sample */
int mm=9;
double xx[20] = { 23,  21,  22,  14,  23,  19,  15,  15,  16 };

/* Second sample */
int nn=14;
double yy[20] = { 27,  27,  22,  21,  27,  28,  19,  20,  31,  24,
   21,  23,  23,  21 };

/* Arrays for the two samples together */
double ww[MAXMN], mrank[MAXMN], mrankx[MAXMN], mranky[MAXMN];
/* Tie-group sizes */
int tgsize[MAXMN];

/* Which sample and which tiegroup */
char type[MAXMN], tgname[MAXMN];


int main(void)
 { int i,j,k,n, r,g,nmn, ntiegroups, nmwdifs,degfree;
   double sum, tiesum,tiecorr, Wval,WExpMean,varnotc,vartc;
   double wmed, zznotc,zztc, pval,pvalwrs, xbar,ybar,zbar;
   int cc1,cc2, m1,m2;  double w1,w2;
   double sx1,sy1, sspool,tt, ttstderr,twid;
   printf ("\nIllustrating nonparametric two-sample methods:\n");
   printf ("(Wilcoxon rank sum and allied procedures - %s.)\n", PROG_DATE);
   printf ("(This is  %s  on the Ma408 Web site.)\n", OUTPUT_NAME);
   /* Display the 1st sample */
   printf ("\nSample I (m=%d):\n", mm);
   for (i=0; i<mm; i++) 
     printf ("  %g", xx[i]);  printf("\n");
   /* Display the 2nd sample */
   printf ("Sample II (n=%d):\n", nn);
   for (i=0; i<nn; i++) 
     printf ("  %g", yy[i]);  printf("\n");

   /* Find joint ranks: First, combine samples together */
   for (r=i=0; i<mm; r++, i++)
    { type[r]='X';  ww[r]=xx[i]; }
   for (j=0; j<nn; r++, j++)
    { type[r]='Y';  ww[r]=yy[j]; }
   nmn=mm+nn;   /* Size of two samples together */
   /* Sort first before finding midranks */
   for (i=0; i<nmn; i++) for (j=i+1; j<nmn; j++)
    { if (ww[i]>ww[j])
       { char cc;  double temp=ww[i];  ww[i]=ww[j];  ww[j]=temp;
         cc=type[i];  type[i]=type[j];  type[j]=cc; }}
   /* Now find midranks */
   for (g=r=0; r<nmn; )   /* r is first value of g-th tie group */
    { double val=ww[r], sumr, midr;
      /* Find the extent of the g-th tiegroup */
      for (i=r+1; i<nmn; i++)  if (ww[i]!=val) break;
      /* The tie group is between r and i-1, inclusive */
      tgsize[g]=n=i-r;  /* Number of elements in tiegroup */
      /* Compute the midranks = average tiegroup rank */
      for (sumr=0, k=r; k<i; k++) sumr+=(k+1);
      midr=(double)sumr/n;
      /* Assign midrank `midr' to the tiegroup */
      /* In ASCII code, 'A'+1='B', 'A'+2='C', etc */
      for (k=r; k<i; k++)
       { mrank[k]=midr;  tgname[k]='A'+g; }
      r=i;  g++;  /* Set for next tiegroup */ }
   ntiegroups=g;  /* Number of tie groups */

   printf ("\nSorted data from both samples:\n");
   printf (" Samp  Value  Tgrp  Ord  Midrank\n");
   for (i=0; i<nmn; i++)
     printf ("  %c    %3g     %c  %3d   %4.1f\n",
       type[i], ww[i], tgname[i], i+1, mrank[i]);
   /* How many non-singleton tie groups? */
   printf ("\n%d tiegroups\n", ntiegroups);
   for (n=g=0; g<ntiegroups; g++) if (tgsize[g]>1) n++;
   printf ("%d tiegroups with more than one tied value\n", n);

   /* Apportion midranks back to the two samples: */
   i=j=0;
   for (k=0; k<nmn; k++)
    { if (type[k]=='X') mrankx[i++]=mrank[k];
      else mranky[j++]=mrank[k]; }
   printf ("Check: End reassortment:  i=%d m=%d   j=%d n=%d\n", i,mm,j,nn);

   printf ("\nThe ranks/midranks in each sample:\n");
   for (sx1=0.0, i=0; i<mm; i++) sx1 += mrankx[i];
   printf ("Sample I (m=%d) (Rank sum=%g):\n", mm,sx1);
   for (i=0; i<mm; i++)
     printf ("  %g", mrankx[i]);  printf("\n");
   for (sy1=0.0, j=0; j<nn; j++) sy1 += mranky[j];
   printf ("Sample II (n=%d) (Rank sum=%g):\n", nn,sy1);
   for (j=0; j<nn; j++)
     printf ("  %g", mranky[j]);  printf("\n");
   printf ("Check: Rsumx+Rsumy=%g:  MN*(MN+1)/2 = %d*%d/2 = %d\n",
     sx1+sy1, nmn,nmn+1, (nmn*(nmn+1))/2);

   /* Compute W statistic based on the second sample */
   Wval=sy1;   /* Rank sum already calculated */
   WExpMean=nn*(nmn+1)/2.0;   /* E(W)  under random permutations */
   printf ("\nObserved Wilcoxon rank-sum statistic based on Sample II:\n"
     "  W=%g   E(W)=%g   W-E(W)=%g\n",  Wval, WExpMean, Wval-WExpMean);
   /* Var(W) WITHOUT tie correction */
   varnotc=(mm*nn*(nmn+1))/12.0;
   zznotc=(Wval-WExpMean)/sqrt(varnotc);
   pval=normpval(zznotc);
   printf ("  Normal approx. W/O tie corr.:  Z=%.3f  P=%.4f (two-sided)\n",
     zznotc,pval);
   /* CORRECT Var(W) WITH tie correction */
   /* Find cubic tiegroup sum for tie correction */
   tiesum=0.0;
   for (g=0; g<ntiegroups; g++)
    { int tg=tgsize[g];  tiesum += (tg-1)*tg*(tg+1); }
   /* Relative tie correction in formula for Var(W) */
   tiecorr=tiesum/((nmn-1)*nmn*(nmn+1));
   printf ("  %d tiegroups   Tiesum: %g   Tiesum correction: %.4f\n",
     ntiegroups, tiesum, tiecorr);
   vartc=varnotc*(1 - tiecorr);     /* True Var(W) with tie correction */
   printf ("  Correct Z-score with tie correction:\n");
   /* Large-sample normal approximation for W and normal P-value */
   zztc=(Wval-WExpMean)/sqrt(vartc);
   pvalwrs=normpval(zztc);
   printf ("  Normal approximation:  Z=%.3f  P=%.4f (two-sided)\n",
     zztc, pvalwrs);

   /* Find Hodges-Lehmann estimator of Y-X: */
   printf ("\nFinding Hodges-Lehmann estimator of Y-X:\n");
   printf ("Collecting and sorting  N = %d*%d = %d  "
     "Mann-Whitney differences:\n", mm,nn,mm*nn);
   /* Collect MW differences into array ww[] */
   for (r=i=0; i<mm; i++) for (j=0; j<nn; j++)
      ww[r++]=yy[j]-xx[i];
   /* Sorted Mann-Whitney differences by `bubble sort' */
   nmwdifs=mm*nn;
   for (i=0; i<nmwdifs; i++) for (j=i+1; j<nmwdifs; j++)
     if (ww[j]<ww[i]) /* If out of order */
      { double tt=ww[i];  ww[i]=ww[j];  ww[j]=tt; }
   /* HL estimator = median value of sorted array[] */
   m1=m2=nmwdifs/2;  if ( (nmwdifs%2)==0) m1--;
   wmed=(ww[m1]+ww[m2])/2;
   /* Find corresponding sample means and sample mean difference: */
   for (sum=0.0, i=0; i<mm; i++) sum+=xx[i];  xbar=sum/mm;
   for (sum=0.0, j=0; j<nn; j++) sum+=yy[j];  ybar=sum/nn;
   zbar=ybar-xbar;
   /* %.2f means round to 2 places beyond the decimal point */
   printf ("HLEst for ``E(Y)-E(X)'' = Median of MW differences\n");
   printf ("  HLEst: %.2f   Ybar-Xbar = %.2f\n",
     wmed, zbar);

   /* Find the HL-associated CI of Y-X: */
   printf ("\nUsing normal approximation to find indices in (1,%d)\n", mm*nn);
   printf ("  for HL-associated 95%% confidence interval:\n");
   /* floor(x) `chops' x, so  floor(2.681)=2 */
   /* ceil(x) `anti-chops' x, so that  ceil(2.681)=3 */
   /* To be conservative, we want to round cc2 up if it is a fraction */
   cc2=(int)ceil((mm*nn)/2.0 + 1.960*sqrt(mm*nn*(nmn+1)/12.0));
   cc1=mm*nn+1-cc2;
   /* Find cc1 and convert to C offsets */
   cc1--;  cc2--;  w1=ww[cc1];  w2=ww[cc2];
   printf ("  C offsets (%d,%d)   95%% CI: (%g, %g):\n", cc1,cc2, w1,w2);
   printf ("The %d sorted Mann-Whitney differences:\n", nmwdifs);
   printf ("(The starred values for median or 95%% CI)\n");
   for (i=0; i<nmwdifs; i++)
    { int starred = (i==cc1 || i==cc2 || i==m1 || i==m2);
      if (i>0 && (i%16)==0)  printf("\n");
      printf ("  %g%s", ww[i], (starred) ? "(*)" : ""); }
   printf ("\n");

   /* For a comparison, find the Student-t 95% confidence interval: */
   /* Find pooled variance estimator */
   sum=0.0;
   for (i=0; i<mm; i++)  sum+=square(xx[i]-xbar);
   for (j=0; j<nn; j++)  sum+=square(yy[j]-ybar);
   sspool=sum/(mm+nn-2);
   ttstderr=sqrt(((1.0/mm) + (1.0/nn))*sspool);
   degfree=mm+nn-2;
   /* Half-width of Student-t 95% CI */
   twid=crit_tdist(degfree,0.05)*ttstderr;
   tt=zbar/ttstderr;
   pval=pvttest(degfree,tt);

   printf ("\nSummary:\n");
   printf ("           Test:Z   P-value   EstDiff    95%% CI\n");
   printf ("WRS:        %5.3f    %6.4f    %.2f    (%.2f, %.2f)\n",
     zztc, pvalwrs, wmed, (double)w1, (double)w2);
   printf ("Student-t:  %5.3f    %6.4f    %.2f    (%.2f, %.2f)   df=%d\n",
     tt, pval, zbar, zbar-twid, zbar+twid, degfree);

   /* No errors to report */
   return 0; }