# exact Monte Carlo simulation of a large sample histogram

(Asked this on MathOverflow, got redirected here)

Say I want to create a histogram of $N$ samples from some simple compactly supported distribution on $\mathbb{R}$, where $N$ is very large, say $N = 10^{30}$. The histogram has $K$ disjoint bins, where $K$ is a more reasonable number like $10$ or $1000$. Obviously it's not feasible for me to directly draw $N$ Monte Carlo samples from my distribution and bin them up to form the histogram. However, it seems to me that there might be some correct method which works by sampling the number of counts in each bin, one at a time, for a total of only $K$ samples. I'm looking for some help finding such a method.

The number of counts in each bin is a binomial random variable with parameters that can easily be calculated by integrating the distribution over the bin interval. So if I have a good way to simulate binomial random variables with large means, I can simulate the number of counts in each bin using only $K$ Monte Carlo draws of a binomial random variable.

The problem is that the total counts in the bins are correlated by the constraint that they must add up to $N$. My method will produce a random number of total counts which will almost certainly not be $N$.

I can think of a couple more sophisticated methods that would avoid this problem, but they create thornier ones - and the bottom line is, I don't know how to prove that any of these heuristic methods are correct.

Can anybody think of an $O(K)$ algorithm that generates a provably correct (or provably nearly correct) histogram sample for this kind of problem? More formally, a correct method for sampling the random vector $H \in \mathbb{Z}^K$ whose entries are the histogram counts? If not, what's the best that I can hope for?

• It seems you can calculate each bin's count's expected value. Why do you need to simulate anything? The simulation will give you a random variable centered around the expected value you can calculate. – Seth May 30 '12 at 20:31
• How small can the bin probabilities get? Will they be within a few orders of magnitude of $10^{30}/K$? – whuber May 30 '12 at 20:33

Simulation (Sheldon Ross) gives an algorithm in section 4.6 by noting that an outcome $X_1,\ldots,X_K$ can be generated in sequence because $X_1$ has a binomial distribution and each conditional distribution $X_i | X_{i-1}, \ldots, X_1$ is also binomial.

Specifically, if the $K$ bin probabilities are $p_1, p_2, \ldots, p_K$, draw $x_1$ from a binomial$(N,p_1)$ distribution, then draw $x_2$ from a binomial$(N-x_1, p_2/(1-p_1))$ distribution, ..., $x_{i+1}$ from a binomial$(N-x_1-\cdots-x_i, p_{i+1}/(1-p_1-\cdots-p_i))$ distribution, etc.

Here is a (recursive) Mathematica implementation.

ClearAll[f];
f[n_, p_] /; Length[p] >= 2 :=
With[{x = RandomInteger[BinomialDistribution[n, First[p]]]},
{x, f[n - x, Rest[p] / (Plus @@ Rest[p])]}];
f[0, p_] := 0 p;
f[n_, p_] /; Length[p] == 1 := n;


Example of use:

x = Block[{p = 1/Range[500], $RecursionLimit = 2000}, f[10^30, p/(Plus @@ p)] // Flatten ];  (4.7 seconds). When the expectations are all greater than 10 or so you can do this much faster by approximating the multinomial as a (degenerate) multivariate normal distribution. • Unless your bins are strangely chosen, the expected count in each one will be$10^{20+}$and so the standard deviations of those counts will be around$10^{10+}$, less than$10^{-10}$of the expectations: do you really need that kind of precision? In a histogram?? – whuber May 30 '12 at 21:08 • Strange as it sounds, this is what I want. I'm trying to do diagnostics on a statistical inference method which appears to be excessively sensitive to small perturbations from the expected values of these histogram bins. I want to feed the algorithm very small (but statistically correct) perturbations and gradually make them larger to see how the performance degrades. – Paul May 31 '12 at 2:12 • Paul, if the histogram bins aren't strangely small themselves, then the normal approximation will be incredibly good. It would take the universe's lifetime of testing to distinguish variates generated by the exact multinomial and constrained normal methods. – whuber May 31 '12 at 16:09 • Thanks, Paul. If you ever find you need better performance, you can use the sequential Binomial algorithm for the low-count bins and then move over to a multivariate normal algorithm for the remaining bins. – whuber Jun 1 '12 at 14:31 • If you're interested in what kind of bizarre application this could be... it's not really a histogram of a distribution on$\mathbb{R}$. The distribution is actually in$\mathbb{R}^3$and it represents a data-cube which records certain properties of ion impacts in a scientific detector. Some parts of the data-cube get lots of counts and others get very few. My interest is in simulating the data-cube realistically. (You don't really see$N= 10^{30}\$ ions, that's just for the debugging process I mentioned above). – Paul Jun 1 '12 at 14:32