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?
 A: 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.
