Testing statistical software

Question

What techniques/approaches are useful in testing statistical software? I'm particularly interested in programs that do parametric estimation using maximum likelihood.

Comparing results to those from other programs or published sources is not always possible since most of the time when I write a program of my own it is because the computation I need is not already implemented in an existing system.

I am not insisting on methods which can guarantee correctness. I would be happy with techniques that can catch some fraction of errors.

Answer 1

One useful technique is monte carlo testing. If there are two algorithms that do the same thing, implement both, feed them random data, and check that (to within a small tolerance for numerical fuzz) they produce the same answer. I've done this several times before:

• I wrote an efficient but hard to implement $O(N\ log\ N)$ implementation of Kendall's Tau B. To test it I wrote a dead-simple 50-line implementation that ran in $O(N^2)$.

• I wrote some code to do ridge regression. The best algorithm for doing this depends on whether you're in the $n > p$ or $p > n$ case, so I needed two algorithms anyhow.

In both of these cases I was implementing relatively well-known techniques in the D programming language (for which no implementation existed), so I also checked a few results against R. Nonetheless, the monte carlo testing caught bugs I never would have caught otherwise.

Another good test is asserts. You may not know exactly what the correct results of your computation should be, but that doesn't mean that you can't perform sanity checks at various stages of the computation. In practice if you have a lot of these in your code and they all pass, then the code is usually right.

Edit: A third method is to feed the algorithm data (synthetic or real) where you know at least approximately what the right answer is, even if you don't know exactly, and see by inspection if the answer is reasonable. For example, you may not know exactly what the estimates of your parameters are, but you may know which ones are supposed to be "big" and which ones are supposed to be "small".

Answer 2

Not sure if this is really an answer to your question, but it is at least tangentially related.

I maintain the Statistics package in Maple. An interesting example of difficult to test code is random sample generation according to different distributions; it is easy to test that no errors are generated, but it is trickier to determine whether the samples that are generated conform to the requested distribution "well enough". Since Maple has both symbolic and numerical features, you can use some of the symbolic features to test the (purely numerical) sample generation:

1. We have implemented a few types of statistical hypothesis testing, one of which is the chi square suitable model test - a chi square test of the numbers of samples in bins determined from the inverse CDF of the given probability distribution. So for example, to test Cauchy distribution sample generation, I run something like

   with(Statistics):
   infolevel[Statistics] := 1:
   distribution := CauchyDistribution(2, 3):
   sample := Sample(distribution, 10^6):
   ChiSquareSuitableModelTest(sample, distribution, 'bins' = 100, 'level' = 0.001);
   

   Because I can generate as large a sample as I like, I can make $\alpha$ pretty small.

2. For distributions with finite moments, I compute on the one hand a number of sample moments, and on the other hand, I symbolically compute the corresponding distribution moments and their standard error. So for e.g. the beta distribution:

   with(Statistics):
   distribution := BetaDistribution(2, 3):
   distributionMoments := Moment~(distribution, [seq(1 .. 10)]);
   standardErrors := StandardError[10^6]~(Moment, distribution, [seq(1..10)]);
   evalf(distributionMoments /~ standardErrors);
   

   This shows a decreasing list of numbers, the last of which is 255.1085766. So for even the 10th moment, the value of the moment is more than 250 times the value of the standard error of the sample moment for a sample of size $10^6$. This means I can implement a test that runs more or less as follows:

   with(Statistics):
   sample := Sample(BetaDistribution(2, 3), 10^6):
   sampleMoments := map2(Moment, sample, [seq(1 .. 10)]);
   distributionMoments := [2/5, 1/5, 4/35, 1/14, 1/21, 1/30, 4/165, 1/55, 2/143, 1/91];
   standardErrors := 
     [1/5000, 1/70000*154^(1/2), 1/210000*894^(1/2), 1/770000*7755^(1/2), 
      1/54600*26^(1/2), 1/210000*266^(1/2), 7/5610000*2771^(1/2), 
      1/1567500*7809^(1/2), 3/5005000*6685^(1/2), 1/9209200*157366^(1/2)];
   deviations := abs~(sampleMoments - distributionMoments) /~ standardErrors;
   

   The numbers in distributionMoments and standardErrors come from the first run above. Now if the sample generation is correct, the numbers in deviations should be relatively small. I assume they are approximately normally distributed (which they aren't really, but it comes close enough - recall these are scaled versions of sample moments, not the samples themselves) and thus I can, for example, flag a case where a deviation is greater than 4 - corresponding to a sample moment that deviates more than four times the standard error from the distribution moment. This is very unlikely to occur at random if the sample generation is good. On the other hand, if the first 10 sample moments match the distribution moments to within less than half a percent, we have a fairly good approximation of the distribution.

The key to why both of these methods work is that the sample generation code and the symbolic code are almost completely disjoint. If there would be overlap between the two, then an error in that overlap could manifest itself both in the sample generation and in its verification, and thus not be caught.

Answer 3

Bruce McCullough had a bit of a cottage industry in assessing statistical software (in the widest sense; he also tested Microsoft Excel. And found it wanting). Two papers that illustrate part of his approach are here and here.

Answer 4

A lot of detail is given by the President of StataCorp, William Gould, in this Stata Journal article.1 It is a very interesting article about quality control of statistical software.