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
p_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
set.seed(1234)
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
var_estimate
> 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.
