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

Question

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)$.

Answer 1

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.$