Exact simulation of a large sample histogram Say I want to create a histogram of $N$ random points 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 points 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 points 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 points 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: Can we assume the number of counts in each bin are large enough that the binomial distribution of counts in the bin will be close to a normal distribution?  And in that case, is it enough to just sample the normal distribution for each bin, add up the samples, and get a total $N'=a\cdot N$ where $a=1\pm O\big({1\over \sqrt K}\big)$ or something like that.  Then just scale all the counts by a factor of $1/a$ to get them to add up to the right thing (scaling the normal distributions keeps them normal).  Am I missing something significant?
A: In short: You want efficient exact simulation from a given multinomial distribution $\text{Mult}(N, p_1,\ldots,p_K)$, where $N$ is very large but $K$ is smallish, for example, $N=10^{30}$ and $K=1000$.
You have already observed that each bin count $C_k$ is binomial with parameters $N$ and $p_k$, where $p_k$ is the probability of a random value from the underlying distribution to be in the $k$th interval.  You have also assumed that the bin probabilities $p_1,\ldots,p_K$ can be calculated ("easily").
What remains is two things.

*

*You correctly observed that the bin counts are correlated.
From your question, it seems that dealing with this is the key
issue you are facing. But this part will be easy and exact.

*How to generate random numbers from binomial
distributions $\text{Bin}(M, p)$ where $M$ is very large.
This part may be easy or hard, depending on your requirements of
speed and accuracy.

You can just generate the bin counts one by one
The individual bin counts are binomial. Their joint distribution is  multinomial with parameters $N$ and $(p_1,\ldots,p_K)$.  So we need to sample a vector of bin counts $\mathbf{C} = (C_1,\ldots,C_K)$ from that distribution. The following algorithm does it.

*

*Initialize $M \leftarrow N$ and $q \leftarrow 1$.

*For $k=1,\ldots,K$, do the following:

*

*Generate the bin count $c_k$ randomly from $\text{Bin}(M, \; p_k/q)$.

*Decrease $M$ by $c_k$ and decrease $q$ by $p_k$.



Why does it work?  Clearly, the first bin count comes from the correct distribution $\text{Bin}(N, p_1)$.  The key observation is that having observed $C_1 = c_1$, the conditional distribution of the remaining bins is multinomial with the remaining count
$N-c_1$ and with the remaining probabilities $p_2,\ldots,p_K$,
scaled to sum to one.
Intuitively: having observed that exactly $c_1$
of your $N$ original points landed on the first interval, the remaining $N-c_1$ points are distributed over the remaining $K-1$ intervals and the relative probabilities of those bins are unchanged.  It is not difficult to prove this, by considering the conditional probabilities of points landing on different bins.
More generally, having observed some of the counts from a multinomial distribution, the conditional distribution of the remaining counts is again multinomial.  (And then the marginal distribution of a single count is binomial.)  See e.g. here in LibreTexts.

The multinomial distribution is also preserved when some of the
counting variables are observed. Specifically, suppose that $(A,B)$ is
a partition of the index set $\{1,2,\ldots,k\}$ into nonempty subsets.
Suppose that $(j_i:i \in B)$ is a sequence of nonnegative integers,
indexed by $B$ such that $j=\sum_{i \in B} j_i \le n$. Let
$p=\sum_{i \in A} p_i$.
The conditional distribution of $(Y_i:i \in A)$
given $(Y_i=j_i:i \in B)$ is multinomial with
parameters $n−j$ and $(p_i/p : i \in A)$.

Generating a random number from a large binomial distribution
How to generate a random number from $\text{Bin}(M, p)$ where $M$ and $Mp$ are large ($M$ being the "remaining" count)? If $M$ were small, you might use the inverse CDF method, like Wikipedia suggests: generate a uniform variate $u$, then add up bin probabilities until you exceed $u$. But this takes $O(M)$ time, so infeasible in your case.
Interestingly, Matlab, Octave and R all fail if we call their binomial random generators with large $M$.  For example Matlab 9.10.0.1602886 (R2021a):
>> binornd(1e5, .2)
ans =
       19909
>> binornd(1e30, .2)
Error using rand
Requested array exceeds the maximum possible variable size.

Given that $M$ is very large, for practical purposes I would just
approximate the binomial with a Gaussian.  If the bin probabilities are in a suitable range, the approximation should be fairly good (but beware of bad approximation when far from the mode).
If that is not an option, I would search the literature for "exact generation of binomial distribution".  At least the article by Kachitvichyanukul and Schmeiser (1989) looks promising; they describe a method they call BTPE (Binomial, Triangle, Parallelogram, Exponential), and run actual experiments up to $N=10^8$. The running time remains quite reasonable. I do not know how well their method scales to $N=10^{30}$.
One more option is to use the inverse CDF method, if you have fast and accurate enough approximation for binomial CDF (or inverse CDF). For example, although R (3.6.3) fails on generating rbinom(1, 1e30, 0.2), it happily computes inverse CDF values, so qbinom(runif(1), 1e30, 0.2) works fine and fast. About the accuracy I do not know, but one could inspect the source.
Reference
Kachitvichyanukul, Voratas; Schmeiser, Bruce W., Algorithm 678: BTPEC: Sampling from the binomial distribution, ACM Trans. Math. Softw. 15, No. 4, 394-397 (1989). ZBL0900.65007.
