Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

1 Answer

up vote 10 down vote accepted

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]] ]

Some notable things about this code:

  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 *)
share|improve this answer
3  
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
1  
@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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.