Sampling Distribution for the Wald-Wolfowitz 2-Sample Runs Test, R

Sampling Distribution for the Wald-Wolfowitz 2-Sample Runs Test

# R version 3.3.2 
# Compute the sampling distribution for the
# two-group Wald-Wolfowitz runs test. 
# (Corrected program, version 2.0)
#
# (Note number of runs can be determined from 
# RLE function, or by visual inspection.)
#
# User sets group sizes here:
n1=50
n2=50
# Define a few functions
partfact = function(n,k) { 
  # compute the partial factorial of n to k
  h=n-k+1
  if(h>n) { return(1) }
  if(h==n) { return(n) }
  if(h<n) { p=prod(seq(h,n)); return(p) }
}
fact = function(k) {
  # compute the factorial of k
  m=1
  if(k>1) { for(i in 2:k) { m=m*i } }
  return(m)
}
# Misc. housekeeping stuff
g=n1+n2          # total number of observations N
p=as.vector(g+1) # define vector of probabilities for successive r values 
sum=0            # clear cumulative prob. accumulator
#
# Compute individual probabilities of r (the number of runs)
# (can't have any less than 2 and any more than N runs)
for(r in 2:g) {     # do for r=1 to N:
    if((r/2) == floor(r/2)) {     # is r even?
        # yes, use r even formula
        k=(r/2)-1
        # compute binomial coefficients
        n=n1-1; d=k; e=partfact(n,d); f=fact(d); a=e/f
        n=n2-1; d=k; e=partfact(n,d); f=fact(d); b=e/f
        n=n1+n2; d=n1; e=partfact(n,d); f=fact(d); c=e/f
        # compute prob. of r
        p[r]=(2*a*b)/c
        }
    else {
        # no, use r odd formula
        k=(r-1)/2
        # compute binomial coeffs.
        n=n1-1; d=k; e=partfact(n,d); f=fact(d); s=e/f
        n=n2-1; d=k-1; e=partfact(n,d); f=fact(d); t=e/f
        n=n2-1; d=k; e=partfact(n,d); f=fact(d); u=e/f
        n=n1-1; d=k-1; e=partfact(n,d); f=fact(d); v=e/f
        n=n1+n2; d=n1; e=partfact(n,d); f=fact(d); w=e/f
        # compute prob. of r
        p[r]=((s*t)+(u*v))/w
        }
    }
# print prob. of r and cumulative prob. distribution
cat("Sampling distribution for Wald-Wolfowitz 2-sample \n")
cat("runs test with n1=", n1, " and n2=", n2, "...\n\n")
for(i in 2:g) {
    sum=sum+p[i]     # compute cum. prob. dist.
    cat(sprintf('For r= %3d: p(r)= %.6e, alpha= %f) \n',i,p[i],sum)) #pretty print results
}
# end.


run \| edit \| history \| help	0