I wanted to ask 1) if I've solved this puzzle problem correctly, and 2) if there is a shorter or more elegant approach.

There are 43 cookies to be given out at random to 10 children.
What is the probability that each child gets at least 2 cookies?

First, notice that this is a classic multinomial setup.

The multinomial distribution can be thought of as giving the probability of observing the outcome $i \in \{ 1, ..., k \}$ coming up $x_i$ times when rolling a $k$-sided die $n$ times, and is a generalization of the Binomial distribution. We have that the probability of each outcome $i$ coming up on a single roll is given by $\pi_i$.

If the so-called die is fair, $\pi_i = \pi_{i'}\; \forall i, i' \in \{ 1, ..., k\}$ or the die is unfair and $\pi_i$ is not necessarily equal to $\pi_i'$. In all cases, we assume $\sum \pi_i = 1$.

Let $X = (X_1, ..., X_k)$ be a $\text{Multinomial}(n, k, \pi)$ distributed where $\pi = (\pi_1, ..., \pi_k)$.

The density of the Multinomial distribution is

\begin{align} P(X = x) & = {n \choose x_1,...,x_k!} \prod_{i=1}^k \pi_i^{x_i} \\\\ & = \frac{n!}{x_1!\cdots x_k!}\pi_1^{x_1} \cdots \pi_k^{x_k}. \end{align}

Here, let $X_1, ..., X_{10}$ denote how many of the $n=43$ cookies each of the $k=10$ children received.

Let $p_n$ denote the probability that all 10 children receive at least 2 cookies each given that $n$ are given out uniformly at random.

Notice that $$ {\small X_1, ..., X_{10} \Bigg\vert \sum_{i=1}^{10} X_i = 43 \sim \text{Multinomial}(n = 43, k = 10, \pi_i = 1/10),} $$ where $\pi_i = 1/10$ indicates that the distribution of cookies is uniformly random (e.g., equal probability) across the 10 children.

The Poisson and Multinomial distributions have an interesting relationship. When the outcomes $X_1, ..., X_k$ are such that $X_i \sim \text{Poisson}(\lambda_i)$, then $\sum_{i=1}^k X_i \sim \text{Poisson}(\sum_{i=1}^k \lambda_i).$

One can easily derive via the definition of conditional probability that $$(X_1, ..., X_k) {\Bigg\vert} \sum_{i=1}^k X_i = N = n \sim \text{Multinomial}(n, k, \pi),$$ where $\pi = \left(\frac{\lambda_1}{\sum \lambda_i}, ..., \frac{\lambda_k}{\sum \lambda_i}\right)$.

\begin{align} P(X = x \Big \vert N = n) & = \frac{P(X = x, N = n)}{P(N = n)} \\\\ & = \frac{P(X = x)}{P(N = n)} \\\\ & = \left( \frac{e^{-\sum \lambda_i} \prod \lambda_i^{x_i}}{\prod x_i!} \right) \Bigg / \left( \frac{e^{-\sum \lambda_i} (\sum \lambda_i)^{n}}{n!} \right) \\\\ & = {n \choose x_1, ..., x_k} \frac{\prod \lambda_i^{x_i}}{\left( \sum \lambda_i \right)^n} \\\\ & = {n \choose x_1, ..., x_k} \left( \frac{\lambda_i}{\sum \lambda_i} \right)^{x_i} \\\\ & \sim \text{Multinomial}(n, k, \pi). \end{align}

Since no information was given to us about how it came about that $n = 43$ cookies were given out, let's assume that it was the result of a $\text{Poisson}(\lambda)$ process. This implies that the $X_i$ are independently and identically distributed as $\text{Poisson}(\lambda/10)$ by a similar argument in the reverse direction.

$$P(X = x \Bigg| \sum X_i = n) = \frac{P(X = x, \sum X_i = n)}{\underbrace{P(\sum X_i = n)}_{\text{Poisson}(\lambda)}},$$ and $X = x \implies \sum X_i = \sum x_i = n$, so $P(X = x, \sum X_i = n) = P(X = x)$ as long as $n = \sum x_i$.

\begin{align} \therefore \;\; P(X = x) & = P(X = x \Bigg| \sum X_i = n) \cdot P(\sum X_i = n) \\\\ & = \left[ {n \choose \pi_1, ..., \pi_k} \prod \pi_i^{x_i} \right] \cdot \left( \frac{e^{-\lambda} \lambda^{\sum x_i}}{\sum x_i! } \right) \\\\ & = \frac{n!}{x_1!\cdots x_k!} \pi_1^{x_1}\cdots \pi_k^{x_k} \left( \frac{e^{-\lambda} \lambda^{n}}{n!} \right). \end{align}

Assume as given in the problem that the cookies are uniformly randomly given out, and $\pi_i = \pi_{i'}\; \forall i, i' \in \{ 1, ..., k \}$; this single probability must be $1/k$ (or in our case, $k = 10$ children).

\begin{align} P(X = x) & = \frac{e^{-\lambda}}{x_1! \cdots x_k!} \left( \frac{\lambda}{k} \right)^{\sum x_i} \\\\ & = \prod_{i=1}^k \frac{e^{-\lambda/k} \left( \frac{\lambda}{k} \right)^{x_i}}{x_i!}, \\\\ \text{a product of the }&\text{density for $k$ iid Poisson}\left(\frac{\lambda}{k}\right) \text{ variables}. \end{align}

We said that if $n$ cookies are given out, then there's a $p_n$ probability that the 10 children all receive 2+ cookies. Then without conditioning on the number of cookies given out, the probability that all kids receive 2+ cookies is given by $$\mathbb{E}\left[\mathbb{E}[p_n | N = n]\right] = \sum_{n=0}^\infty \frac{e^{-\lambda} \lambda^n}{n!} p_n.$$

On the other hand, since the $X_1, ..., X_k$ are iid Poisson, the probability that all 10 kids receive 2+ cookies is $P(X_1 \geq 2)^{10} = (1-P(X_1 \leq 1))^{10} = (1-e^{-\frac{\lambda}{10}} - \frac{\lambda}{10}e^{-\frac{\lambda}{10}})^{10}$.

Now we have that $$ \sum_{n=0}^\infty \frac{e^{-\lambda} \lambda^n}{n!} p_n = \left[1-e^{-\frac{\lambda}{10}} - \frac{\lambda}{10}e^{-\frac{\lambda}{10}}\right]^{10}. $$

Now, recall that $e^{x}$ has a series expansion. If we multiply both sides by $e^{\lambda}$ and by $43!$, we should find that the coefficient on the left-hand-side series for $\lambda^{43}$ is just $p_n$, and the right-hand-side series coefficient for $\lambda^{43}$ gives an expression we can evaluate for $p_{43}$. They must be equal, because in order for two convergent power series to coincide on a non-empty interval, their coefficients must be equal.

So what we're going to do is expand the function $e^{x} = \sum_{t=0}^\infty \frac{x^t}{t!}$ wherever we see it on the modified right-hand-side and use that to identify the coefficient of $\lambda^{43}$, which is $p_{43}$.

However, this is a bit hard to do by hand, so we'll use the symbolic algebra library sympy in Python to do it for us.

# Python code
from sympy import symbols, exp, factorial, series, Pow

lambda_ = symbols('lambda')

inner_expression = 1 - exp(-lambda_/10) - (lambda_/10)*exp(-lambda_/10)

raised_expression = inner_expression**10 

complete_expression = exp(lambda_) * raised_expression

expanded_series = series(complete_expression, x=lambda_, n=44).removeO()

coeff_lambda_43 = expanded_series.coeff(lambda_**43)

p_43 = factorial(43) * coeff_lambda_43

$$\frac{38360235213946776318553037176114920309}{78125000000000000000000000000000000000} \approx 0.491 $$

Thus we conclude that if 43 cookies are given out to 10 children uniformly at random, then the probability that each child receives at least 2 cookies is $\approx .491.$

Let's see if we can confirm that via a simple simulation in R:

# R code

num_trials <- 100000  # Number of simulations
num_children <- 10    # Number of children
num_cookies <- 43     # Number of cookies

results <- replicate(num_trials, {
  cookies <- sample(1:num_children, num_cookies, replace = TRUE)
  counts <- table(factor(cookies, levels = 1:10))
  all(counts >= 2)

prob_estimate <- mean(results)
var_estimate <- var(results) / num_trials

> prob_estimate
[1] 0.49178
> var_estimate
[1] 2.499349e-06

The point-estimate is off by 0.0007689893. I'm not sure if the way I've calculated the variance is appropriate.


3 Answers 3


This idea is legitimate. We can summarize the series calculation as finding the coefficient of $\lambda^{43}/43!$ in $$F(\lambda) = e^{\lambda}\left(\sum_{k\ge 2} e^{-\lambda/10} \frac{(\lambda/10)^k}{k!}\right)^{10}$$ We can simplify a little bit to get $$F(\lambda) = \left(\sum_{k\ge 2} \frac{(\lambda/10)^k}{k!}\right)^{10}$$ and from here, we can explain why this works directly. To get a term with $\lambda^{43}$ from this tenth power, we must take $\lambda^{k_i}$ from the $i^{\text{th}}$ factor for $i=1,\dots,10$ where $k_i \ge 2$ for all $i$ and $k_1 + k_2 + \dots + k_{10} = 43$, which corresponds to the case where we distribute $k_i$ pieces of candy to the $i^{\text{th}}$ child. Such a term is more precisely equal to $\frac{(\lambda/10)^{43}}{k_1!\,k_2! \cdots k_{10}!}$, so the coefficient of $\lambda^{43}/43!$ in that term is exactly equal to $\frac{43!}{k_1!\,k_2!\cdots k_{10}!}(\frac1{10})^{43}$, the multinomial probability of this case.

(From this point of view, we don't have to think of $\lambda$ as the rate of a Poisson random variable at all - we can simply think of this as an exponential generating function for the probability we wanted.)

As observed in the question, the probability that a Poisson is at least $2$ is $1$ minus the probability that it's $0$ or $1$, which lets us simplify the expression for $F(\lambda)$ somewhat: we can write $F(\lambda) = \left(e^{\lambda/10} - 1 - \lambda/10\right)^{10}$.

Another way to find the coefficient of $\lambda^{43}/43!$ in $F(\lambda)$ is to take the $43^{\text{rd}}$ derivative with respect to $\lambda$, and then set $\lambda=0$. This is still not something we'd want to do by hand, but it's easier to enter into most computer algebra systems. For example, both Mathematica and Wolfram Alpha accept the following:

D[(Exp[λ/10] - 1 - λ/10)^10, {λ,43}] /. λ -> 0

Another approach is the inclusion-exclusion one: this gives a probability of $$\sum_{a=0}^{10} \sum_{b=0}^{10-a} (-1)^{a+b} \binom{10}{a,b,10-a-b} \frac{43!}{(43-b)!} \left(\frac1{10}\right)^b \left(1-\frac{a+b}{10}\right)^{43-b}$$ where the idea is that the $(a,b)$ term accounts for cases where we pick $a$ children to get $0$ cookies and $b$ children to get $1$ cookie (and leave the remaining children undetermined). So we must decide which children are playing which role in $\binom{10}{a,b,10-a-b}$ ways, which cookies are going to the $b$ children with just one cookie in $\frac{43!}{(43-b)!}$ ways, then multiply by the probability that those cookies are doing what we said they would, $\left(\frac1{10}\right)^b$, and the probability that no other cookies go to these $a+b$ children, $\left(1-\frac{a+b}{10}\right)^{43-b}$.


Randomly distributing $m$ cookies to $n$ children, one cookie at a time, let $f(a,b,r)$ be the probability that each child ends up with at least two cookies, assuming

  • There are $r$ cookies remaining to be distributed.$\\[4pt]$
  • There are $a$ children who have so far received no cookies.$\\[4pt]$
  • There are $b$ children who have so far received exactly one cookie.$\\[4pt]$
  • There are $n-a-b$ children who have so far received at least two cookies.$\\[4pt]$

For $m=43,n=10$, we want to find $f(10,0,43)$.

Shown below is a recursive solution, implemented in Maple.

enter image description here

Equivalent Python code:

import sympy as sp
from functools import lru_cache

n = 10

def f(a, b, r):
    if a == 0 and b == 0:
        return 1
    elif r < 2 * a + b:
        return 0
    elif a == 0:
        return (b / n) * f(0, b - 1, r - 1) + ((n - b) / n) * f(a, b, r - 1)
    elif b == 0:
        return (a / n) * f(a - 1, 1, r - 1) + ((n - a) / n) * f(a, b, r - 1)
        return (a / n) * f(a - 1, b + 1, r - 1) + (b / n) * f(a, b - 1, r - 1) + ((n - a - b) / n) * f(a, b, r - 1)

result = f(10, 0, 43)
  • $\begingroup$ Confirmed: I ran the program of zenodo.org/records/10542669 for distributing 43 balls into 10 urns. The probability for the "most even" case where 7 children get each 4 and 3 children get each 5 cookies is 9.14756463421271E-5 . The probability that no child gets no cookie and no child gets 1 cookie is 38360235213946776318553037176114920309/78125000000000000000000000000000000000 = 0.4910110107385187 ... $\endgroup$ Commented May 13 at 9:22
  • 3
    $\begingroup$ This is a good solution. You should write out the recurrence mathematically because the maple code is harder to read. $\endgroup$
    – qwr
    Commented May 13 at 14:30

General case:

Having $m$ cookies, $n$ children, and minimum number of $r$ cookies for each child, the probability $P_r(m,n)$ is given by

$$\color{blue}{P_r(m,n)=\frac{n!S_r(m,n)}{n^m}} \tag{1}$$

where $S_r$ denotes the $r$-associated Stirling number of the second kind [1], which can be efficiently computed using the following recursive formula:

$$S_r(m+1, n)=n\ S_r(m, n)+\binom{m}{r-1}S_r(m-r+1, n-1) $$

Your problem

In your problem, $m=43$,$n=10$, and $r=2$. The 2-associated numbers are studied here OEIS|A008299, with the recursion $$S_2(m+1, n)=n\ S_2(m, n)+mS_2(m-1, n-1),$$

which yields using (1):

$$\color{blue}{P_2(m+1, n)=P_2(m, n)+\frac{m}{n} \left ( 1-\frac{1}{n} \right )^{m-1}P_2(m-1, n-1)}$$

with boundary conditions $$P_2(m,n)=0, m<2n, P_2(m,1)=1, m\ge2$$ $$P_2(0,1)=P_2(1,1)=0.$$ The python code is given below

import numpy
n_max, m_max = 10, 43
p = numpy.zeros((m_max+1, n_max+1))
p[2:,1] = 1
for n in range(2,n_max+1):
    for m in range((2*n)-1, m_max):
        p[m+1,n] = p[m, n] + (m/n) * p[m-1, n-1] * (1-1/n)**(m-1)

print(p[m_max, n_max])

which returns $$\color{red}{0.49101101073851994}.$$

  • $\begingroup$ Please don't post code as pictures. $\endgroup$
    – qwr
    Commented May 13 at 14:30
  • 1
    $\begingroup$ @qwr I just provided the code. $\endgroup$
    – Amir
    Commented May 13 at 15:47

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .