# How to sample $n$ observations from a multinomial distribution using binomial (or poisson) sampling?

## Context

I have $$n$$ observations which I'd like to sample with replacement for the purpose of bootstrap. A way to think about it is that we have a multinomial distribution with $$n$$ classes and that we'd like to draw from it $$n$$ items.

However, due to computational context, it's easier for me to use either a binomial or poisson distribution. So an alternative (known) way to do approximate sampling will be to count how many times we wish to have each item, but drawing a Binomial($$n$$, $$\frac{1}{n}$$) frequency for each of the items. When $$n$$ is large, we can approximate this using Poisson(1) (i.e.: poisson bootstrap).

## The question

I'd like to use these sampling to get exact multinomial distribution, instead of an approximate one. How can this be achieved using either binomial or poisson distribution?

The issue is that since we use, for each of the items, Binomial($$n$$, $$\frac{1}{n}$$), we have a non-zero chance that the sum of these binoms will be different than $$n$$. If it is smaller then n, then we can just use the values up to the $$n-1$$ item, and for the $$n$$th item just use $$n$$ minus the sum of the previous draws. But if the sum is larger than $$n$$, this method won't work.

An idea I had here is to either use some sort of rejection sampling, by which whenever the sample is not exactly $$n$$, to just try it again (which feels VERY expensive).

Another method I thought about was to randomly shuffle the order of the $$n$$ items, and then take the first $$k$$ items, so that their sum would be equal to $$n$$, and if the $$n$$ items would sum to a smaller number than $$n$$, then to pick the $$n$$th value and have it be n minus the sum of the values up to that point.

The issue is that I don't have a proof that this method would yield a valid multinomial distribution (I suspect it will, but no proof).

Any suggestions/references would be lovely. Thanks upfront!

• Are you seeking a method for general multinomial, with probabilities $p_1, p_2, ..., p_k$ which sum to $1$ or a uniform multinomial ($p_1=p_2=...=p_k = 1/k$)? It sounds like the latter, in which case you should really make that explicit. Sep 5, 2021 at 1:22
• I have an answer that uses only binomial calls, either way but it would be simpler to know which case to answer. Sep 5, 2021 at 3:05
• Nevermind, I did both. Sep 5, 2021 at 3:14

You can do it by progressing conditionally through the categories. I'm going to work from the last category backward (for a particular reason) but it can be done in any order as long as you're consistent in how you go about it.

Equal probability case:

Sample the count in the $$n$$-th category, $$X_n\sim\text{bin}(n,\frac{1}{n})$$.

Sample the count in the $$n-1$$-th category conditionally on $$X_n=x_n$$, i.e. $$X_{n-1}\sim\text{bin}(n-x_n,\frac{1}{n-1})$$.

Sample the count in the $$n-2$$-th category conditionally on $$X_{n-1}+X_n=s_{n-1:n}$$, i.e. $$X_{n-2}\sim\text{bin}(n-s_{n-1:n},\frac{1}{n-2})$$, ... and so on. Naturally, if $$s_{i:n}$$ hits $$n$$, all lower $$X$$'s are $$0$$ for that sample.

Unequal probability case:

Sample the count in the $$n$$-th category, $$X_n\sim\text{bin}(n,p_n)$$.

Sample the count in the $$n-1$$-th category conditionally on $$X_n=x_n$$, i.e. compute $$p_{n-1:n} = \frac{p_{n-1}}{{1-p_n}}$$ (i.e. scaling up the remaining probabilities since we no longer have the last category) and draw $$X_{n-1}\sim\text{bin}(n-x_n,p_{n-1:n})$$.

Sample the count in the $$n-2$$-th category conditionally on $$X_{n-1}+X_n=s_{n-1:n}$$, i.e. compute $$p_{n-2:n} = \frac{p_{n-2}}{{1-p_{n-1}-p_n}}$$ and draw $$X_{n-2}\sim\text{bin}(n-s_{n-1:n},p_{n-2:n})$$, ... and so on. Naturally, if $$s_{i:n}$$ hits $$n$$, all lower $$X$$'s are $$0$$ for that sample.

So you progress through, adjusting both the "n" and "p" to account for the conditioning on what categories are already drawn and what's left to draw from.

If there are a very large number of categories and so typically very small values of $$p_i$$ you may need to pay attention to numerical error in the scaling of those $$p$$ values as you progress.

• Thanks @glen_b, it sounds like it makes sense. Is there a reference for this process? (Maybe with a proof?!) Sep 5, 2021 at 15:41
• I don't have a reference, no -- but you might find it among the exercises in an undergrad mathematical statistics textbook. Can't do it now but I'll try to come back later and outline a proof. It is, presumably clear that in respect of one category ("in category A" vs "not in category A") the distribution is binomial. Sep 5, 2021 at 23:10
• Heh, Just decided to try a search in a spare moment and there was something on it in the first hit. The algorithm relies directly on this fact: online.stat.psu.edu/stat504/lesson/1/1.7/1.7.6 (excluding the last paragraph) ... is that sufficient for your purpose? The second hit was this, on which see the subsection on "Conditioning". Sep 5, 2021 at 23:22
• As expected, the first undergrad textbook that I checked (Mendenhall et al) left it entirely as an exercise; this was why I was suggesting a reference might be difficult. Sep 6, 2021 at 0:47
• This is very cool, thanks @Glen_b ! I wonder how much variation there could be about this. For example, if you can only do the sampling upfront for each of the categories, but then wanted to bing the conditionality in somehow, I wonder how much that could be done. For example, if the number of categories is very large, it's clear that the negative correlation between the items becomes very small, so one could sample independently each item. But then the question is if there is some correction to do to the sum of results so to bring back at least some of the dependance. Sep 6, 2021 at 8:22