# What is the statistical distribution of the number of times a number i would appear after M trials of a poisson distribution?

I want to generate $$M$$ random numbers following a poisson distribution with mean $$k$$, and then write down how many times the number 0 appears, how many times the number 1 appeared, and so on, until some number $$n$$, which is the higher number I'm interested in.

Let's call $$C(j)$$ to the number of times the number $$j$$ would be generated after generating $$M$$ random numbers.

What is the probability distribution of $$C(j)$$?

The idea is, instead of generating $$M$$ random numbers (because $$M$$ could be very high), I'd like to generate instead $$n+1$$ random numbers that would match valid numbers of times each number from $$0$$ to $$n$$ would appear after generating $$M$$ numbers following $$P(k)$$.

You ask about the order statistics of any finite discrete probability distribution. Remembering that the order statistics form a Markov chain, we can immediately code a two-line solution.

Let the possible values of your distribution be $$x_0 \lt x_1 \lt \cdots \lt x_{n}$$ with corresponding probabilities $$p_i\gt 0,$$ $$i = 0, 1, 2, \ldots, n.$$ In a sample of size $$M,$$ let $$M_i$$ ($$C(i)$$ in the question) be the number of times $$x_i$$ appears: its multiplicity.

1. The multiplicity of $$x_0,$$ $$M_0,$$ is Binomial$$(M, p_0).$$

2. Conditional on $$M_0 + M_1 + \cdots + M_k = d,$$ with $$0 \le d \le M,$$ the multiplicity $$M_{k+1}$$ also is Binomial. Its parameters must therefore be $$M-d$$ (the number of remaining values in the sample with values $$k+1$$ or larger) and $$p_{k+1} / P_{k+1}$$ (the relative chance that any single value in this set equals $$x_{k+1}$$), where $$P_{k+1}$$ is the value of the survival function at $$x_k:$$ $$P_{k+1} = 1 - (p_0 + p_1 + \cdots + p_k) = p_n + p_{n-1} + \cdots + p_{k+1}.$$

The latter expression is the right way to compute $$P_{k+1}$$ to avoid double-precision errors: start from $$p_n,$$ add $$p_{n-1},$$ and so forth. (Otherwise, accumulated errors can cause $$p_{k+1}/P_{k+1}$$ to exceed $$1,$$ which will choke most Binomial random number generators.)

Here is an R implementation. It returns the vector $$(M_0, M_1, \ldots, M_n)$$ when given $$M$$ and the vector $$p = (p_0, p_1, \ldots, p_n).$$

s <- function(m, p) {
# Initialization.
P <- rev(cumsum(rev(p)))
k <- rep(0, length(p))

# The algorithm.
for (i in seq_along(p)) {
k[i] <- rbinom(1, m, p[i]/P[i])
m <- m - k[i]
}
k
}


For instance, initializing p to the values of a Poisson$$(4)$$ distribution at $$x = 0, 1, \ldots, 102,$$ and tacking on the remaining probability (about $$8\times 10^{-106})$$ to represent any values greater than $$102,$$ and setting M to 1e100 (a sample of one googol values) produces this result:

m <- floor(lambda + 50 * sqrt(lambda)) # Way, way out in the tail...
p <- c(dpois(seq(0, m - 1), lambda), ppois(m - 1, lambda, lower.tail = FALSE))
set.seed(17)                           # For reproducibility
s(1e100, p)

 [1] 1.831564e+98 7.326256e+98 1.465251e+99 1.953668e+99 1.953668e+99 1.562935e+99 1.041956e+99
[8] 5.954036e+98 2.977018e+98 1.323119e+98 5.292477e+97 1.924537e+97 6.415123e+96
...
[92] 8.304320e+12 3.610576e+11 1.552958e+10 6.608275e+08 2.782288e+07 1.160200e+06 4.746200e+04
[99] 1.989000e+03 9.100000e+01 2.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00


The actual order of the values $$x_i$$ doesn't even matter. (The code shows they are not used in the algorithm.) But for maximal accuracy, index your values $$x_i$$ so that $$p_n \le p_{n-1} \le \cdots \le p_0:$$ that is, place the low-probability outcomes last. This will maintain the greatest possible precision in the vector $$P.$$