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


1 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]

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.