How to actually compute the statistical distance?

I have an algorithm for discrete sampling that I've defined. This is to say that $$A$$, upon input some pdf $$p(x)$$, outputs samples from the pdf $$p(x)$$. There are some nuances (you have to make the support of $$p$$ finite, so you're really only sampling from some distribution within some statistical distance $$SD$$ of the distribution of $$p$$), but nothing too complicated. By statistical distance I imean the "Total Variation Distance":

$$d(D_1,D_2) = \sum_{x\in\mathsf{supp}(D_1)\cup\mathsf{supp}(D_2)}\lvert \Pr[D_1 =x] - \Pr[D_2 = x]\rvert$$

Imagine I can (theoretically) prove that this process truly does sample from a distribution that's within $$SD$$ of the true distribution, and I want to implement this sampler. What are the typical ways to show that $$A$$ produces the distribution that it theoretically should? Given a perfect implementation, this wouldn't be an issue, but I'm hoping to use this to catch bugs.

There are easy naive things to do --- for example, draw some number of samples $$(x_1,\dots,x_n)$$, treat this as an empirical distribution, and compute the statistical distance between that and the (true) pdf $$p(x)$$. My worry about this is that the Statistical Distance problem is complete for the complexity class Statistical Zero Knowledge (see chapter 4 of Vadhan's thesis for a survey), which might hint that this is a poor method of finding bugs.

Of course:

1. It could be that my instance of the problem is easy

2. My problem is presented in a slightly different form (SZK is a circuit class I believe)

Still, this makes me worry that checking the statistical distance of $$A$$ from $$p(x)$$ is an intractable problem. Is this true? Are there better algorithms for it than the naive one I mentioned above?

• please explicitly define statistical distance to make your question self contained – kodlu May 23 '19 at 22:24
• As evidence that this question may be a trifle too broad for this site, I present Modern Challenges in Distribution Testing, submitted in partial fulfillment of the requirements for the degree of PhD in EECS at MIT just nine months ago. – Squeamish Ossifrage May 24 '19 at 3:50
• Indeed @SqueamishOssifrage is right. There was earlier work by Tugkan Batu on this as well. maths.lse.ac.uk/Personal/batu/respub.html#pub – kodlu May 24 '19 at 4:34
• @SqueamishOssifrage glancing through that, it appears that testing if a distribution is $\epsilon$-close to a known distribution (with both distributions over $[n]$ is (very roughly) $\mathsf{poly}(n, \epsilon^{-1})$. Since I'm interested in $\epsilon$ being cryptographicly small (roughly $2^{-60}$), this makes it seems like its computationally infeasible. – Mark May 24 '19 at 5:19

2. Factor your code into a nondeterministic uniform bit sampler and deterministic logic.

• You can use the full suite of standard software testing techniques for to test deterministic part.
• You can use automated formal verification tools to prove the deterministic part is correctly implemented.
• You can use well-studied cryptography—also with deterministic known-answer test vectors—seeded by observations from well-studied unpredictable physical models for the uniform bit sampler.

‘OK,’ you say, ‘but Shirley, you should be able to also use empirical hypothesis tests of random samples—after all, the goal is to test that a random sampler for a distribution is correct, so why not run it and apply a standard empirical hypothesis test of distribution to the samples? Specifically, we can consider the null hypothesis that the software correctly samples from the intended distribution, and apply a hypothesis test with a specified statistical significance level for that null hypothesis!’

This is not wrong. But there are several issues with it.

1. Empirical hypothesis tests of distributions are qualitatively different from the normal kinds of automatic tests you might be familiar with in software engineering, because they necessarily have false alarms:

• A test for deterministic software properties rejects certain samples, because correct software cannot yield them—thus it is a bug if you encounter them.

• A random test for the distribution of a null hypothesis rejects certain samples, even though correct software can yield them—the conceit, of course, is that buggy software will yield these rejected samples with much higher probability than correct software.

What do you do, then, when your test suite rejects the null hypothesis and raises an alarm? Well, it might be a bug, or it might just be a fluke! Better run it repeatedly to make sure that the rate of alarms is exactly what the significance level was set to—so maybe we should run a hypothesis test for the null hypothesis of iid Bernoulli trials weighted by the significance level, but what do you do if that test fails? Well, it might be a bug, or it might just be a fluke…

2. For these empirical tests to be useful, they need high statistical power to detect the bugs that might be in your program. So you need to hypothesize plausible bugs, and confirm that the tests have high statistical power to detect those bugs.

• You could hypothesize specific bugs that occurred during development, and relatively easily make tests for them.

But, you could also just test those in the deterministic part of your program.

• You could hypothesize bounds on statistical distances: null hypothesis, $$D(P \mathbin\| Q) \leq \varepsilon$$; alternative hypothesis, $$D(P \mathbin\| Q) \geq \varepsilon + \delta$$. Then you can choose hypothesis tests that are guaranteed to have certain false alarm rates and statistical power for those alternative hypotheses. You could even choose the two bounds to use different statistical distances, e.g. the null hypothesis might be bounded KL divergence while the alternative hypothesis is a bug of excess Hellinger distance.

But, the numbers of samples required for these tests to attain prescribed false alarm rates and statistical power are usually superlinear in $$1/\varepsilon$$ and $$1/\delta$$—if they exist at all with such guarantees, which for some statistical distances $$D$$ they do not!

• You could just apply standard hypothesis tests and blindly hope they'll catch bugs. Actually in practice this works quite well during development to catch obvious bugs early on.

But, it doesn't get you very far.

3. As your system gets bigger, you will be tempted to expand the tests you do.

• With deterministic property tests, adding a test case to a test suite just means each run of the test suite takes a little more time.

• With $$n$$ random hypothesis tests, if each test case has a false alarm rate $$\alpha$$, the test suite has a false alarm rate $$1 - (1 - \alpha)^n$$ which rapidly approaches $$1$$i.e., false alarms nearly all the time—as $$n$$ increases. What can you do when you add a new test?

• Leave the false alarm rate of the test cases alone.

But, then if you had calibrated the test cases so that the test suite has a prescribed false alarm rate, that rate will go up as you add more tests—even if the software is working perfectly!

• Whenever you add a test case, recalibate all the test cases to have false alarm rate $$\alpha' < \alpha$$ so that $$1 - (1 - \alpha)^n = 1 - (1 - \alpha')^{n + 1}$$, i.e. choose $$\alpha' = 1 - (1 - \alpha)^{n/(n + 1)}$$.

But, the number of samples needed in each test case to maintain the same statistical power at lower false alarm rates tends to grow rapidly with $$n$$.

• Give up on a notion of ‘tests passing’, and instead record a detailed history of test case runs on which you can do meta-analyses later.

But, software people don't like it when the answer from a test suite isn't just PASSED or FAILED, and none of the standard tools like Jenkins are designed to cope with such nonbinary fluidity of notions of test results, as is typical for the white cismale culture in Silicon Valley that dominates software engineering tooling.

4. My name is Squeamish, not Shirley.

5. If we're talking about cryptography, the statistical distances are generally going to be so small you don't have a snowball's chance in a Manhattan trash can on a summer day of detecting them anyway. So go back to step 1: prove your method correct.

Of course, all of this depends on more details about what you're doing. Maybe, for your application and your distribution and your sampler and your type of statistical distance and your bound $$\varepsilon$$ on magnitude of statistical distance and your tolerance $$\delta$$ on magnitude of statistical distance and your development process and your toolchain, you can actually address all these concerns, in which case, great! Go for it!

Use Pearson's chi-squared test. (The degrees of freedom being one less than the number of possible outcomes.)

Statistical distance metrics will give you a number, but which metric should you use? How are you supposed to interpret distance defined using that metric? How large of a distance is too much?

A $$\chi^2$$ test is basically a way to interpret distance (a weighted sum of squared errors) in terms of a p-value. In this case, you can interpret that p-value to be the probability that a correctly implemented sampling algorithm would produce output like that which you observe. (With at least as much "weighted sum of square errors" as observed.)

You choose how many samples you will be using. Then calculate the number of times you expect to see each individual outcome. You sum together each $$(observed - expected)^2 / expected$$ term. The resulting sum should follow a Chi-Squared distribution as long as each $$expected$$ is sufficiently large.

If the expected number of occurrences for any given outcome is too small, either merge low-probability outcomes (pooling) or increase the number of samples taken. Different rules of thumb exist for the minimum number of expected occurences (5, 10, 20, 100). Since you can always algorithmically generate more samples, I suggest 100-1000. (Some sources recommend low minimums because it's normally much harder to get as many samples as you want and because the chi-squared distribution approximation is still good enough.))

With at least 1000 expected occurrences per outcome, the test should work for any arbitrary distribution.

In the context of cryptography, however, empirical tests should not be relied upon. We are concerned with very tiny biases. To get that level of precision, you would need a sample size too large to actually generate.

More importantly, they can't be relied upon because no test can assure you that an algorithm will be safe under adversarial conditions.

Failing a statistical tests should be interpreted as a red flag. Passing such tests, however, is not assurance that an algrithm is safe. You can only detect extreme errors, such as implementation issues. Everything else is indistinguishable from "not extremely broken".

Many implementation errors, however, will not result in statistical tests failing. We're dealing with random numbers. Serious flaws can slip by humans and machines because output will still look random. Instead you should focus 99.99% of your energy into analyzing source code.

Fortunately, the type of algorithm you described isn't hard to implement. I can think of six "well known" algorithms that will do the job. Each requires only the ability to sample from a uniform distribution. The most complicated of these algorithms, the alias method, is still fairly simple.

• While algorithms like the alias method (and other generic sampling methods) aren't that difficult on their own, I have a certain extension of them that I'm trying to evaluate the cost of. The theoretical construction is sound, but an asymptotic analysis of the cost is both less interesting in practice, and difficult to make tight due to some details of the construction. For this reason I'm also leaning on an implementation, and am trying to figure out when I'll know my implementation is "right" besides when major bugs are gone. Of course, this implementation is solely for research purposes. – Mark May 24 '19 at 5:22
• Suppose $X_i \sim P$ are iid samples, and you compute $\chi^2 = n\sum_x (\#\{i : X_i = x\}/n - Q(x))^2/Q(x)$, where $D(P \mathbin\| Q) < \varepsilon$ but $P \ne Q$. What's the distribution of $\chi^2$ (say, the limit as the number of samples grows without bound)? The answer is easy when $P = Q$, as you stated, but how does this help when $P$ and $Q$ are merely close with a quantified bound, yet not the same? – Squeamish Ossifrage May 24 '19 at 7:07