# Wald–Wolfowitz Runs Test

Does Mathematica 8 implement Wald–Wolfowitz runs test for randomness? I can't find it in the documentation. I would like to test some fit residuals.

-
– cormullion May 12 '12 at 9:05

This is not built into the system as far as I can tell. However, it's a really good example of how to construct a fairly sophisticated function using just a few lines of code. I've taken the definition directly from the Wikipedia page, and applied it to Sign[data] on the assumption that if you are testing residuals from a regression, you are interested in runs of positive and negative values.

runstest[x:{__?NumericQ}] :=
Module[{n = Length[x], res = Sign[x], n1, n2, runs, mu, sig2, y},
{n1, n2} = Last /@ Tally[res];
runs = Split[res];
mu = 2 n1 n2/n + 1;
sig2 = (mu - 1) (mu - 2)/(n - 1);
NProbability[y <= Length[runs],
y \[Distributed] NormalDistribution[N@mu, N@sig2]] ]


1. The pattern {__?NumericQ}. This ensures that the function only works on vectors of numeric data.
2. The use of local variables via Module, some of which can be set initially in the first argument of Module, while others can wait till later, perhaps because they depend on other local variables.
3. The use of Tally and Split, built-in functions that automatically divide lists into sublists of like elements.
4. The use of NProbability to work out the one-tailed test. Alternatively the last line could be something like N@CDF[NormalDistribution[mu, sig2], Length[runs]]: it gives the same answer except where the p-value is so tiny that numerical error creeps in (but seriously, who cares if it's $10^{-16}$ or $10^{-20}$?). It's possible that you want to express the test a different way, in which case that line is easy to modify.

Here is some test data that should definitely be rejected as random because of strong serial correlation:

testdata1 =
FoldList[0.98 #1 + #2 &, RandomReal[],
RandomVariate[NormalDistribution[0, 0.5], 100]];


And here is some data that should probably pass the test.

testdata2 = RandomVariate[NormalDistribution[0, 0.5], 100];


Sure enough:

runstest[testdata1]

(* -7.90479*10^-14 *)

runstest[testdata2]

(* 0.451757 *)

-
There is the nice identity NProbability[y <= Length[runs], y \[Distributed] NormalDistribution[N@mu, N@sig2]] == 1 - NProbability[y > Length[runs], y \[Distributed] NormalDistribution[N@mu, N@sig2]] you could use to simplify the last line of your function; it is also much more numerically sound to avoid the subtraction of near-equal quantities whenever possible. – J. M. May 14 '12 at 11:07
@Verbeia: Why do you reject the two testdata? What is the expected value for an accettable random dataset? Thanks! – meriens Apr 5 '13 at 8:17
@meriens you have twice now posted answers that are really comments. I converted this one for you. If you will participate just a little more you will soon have the comment everywhere privilege. – Mr.Wizard Apr 5 '13 at 8:41