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.
