# M buckets out of N have at least one ball

We have N buckets and we start filling them randomly with balls. At the end we know that we have exactly M buckets that have at least one ball in them. Both N and M are given and M < N.

What is the CDF of the number of balls that were used? Or at least the mean and variance?

• Can N equal M? Say we use N balls and each bucket gets a ball – bjschoenfeld Aug 30 '19 at 2:35
• If N=M, then I think there is no solution. Might as well be infinite number of balls. – Denisevi4 Aug 30 '19 at 2:40
• I think there is something strange about the question. Even with M < N, you could still have many more balls than N, just one of the buckets happenend not to get a ball. If you had K balls, I could tell you the pmf over M. Given M, you could then choose the K that gave the highest probabilty for M, but I don't know if there is a distribution over K... – bjschoenfeld Aug 30 '19 at 2:57
• CDF distribution is conditional on the observation M. Say M=99 while N=100. If number of balls is K. Obviously K >= 99. If K == 99 we should be able to calculate probability of such event that each of the 99 buckets is filled exactly once. I'm sure it's tiny. If we increase K , probability of M observation would increase for a while. But if K == infinity, probability is zero. So, there must be a maximum for K and must have some distribution. – Denisevi4 Aug 30 '19 at 3:14
• There is no CDF, unless you are interested in the CDF as a function of the number of balls. One can compute the likelihood of $M.$ – whuber Aug 30 '19 at 12:10

Your question is a bit ambiguous, but I'm going to assume that you keep filling the buckets until you get the desired number of balls, and then you count the number of balls you used (i.e., I'm going to interpret this as a "stopping rule"). Under this interpretation, your question is a generalisation of the coupon-collector problem, and it can be solved using the classical occupancy distribution. Instead of using your notation, I'm going to use standard notation for this problem, where we have $$m$$ bins (buckets) and we allocate balls until there are $$m_*$$ occupied buckets.

As a preliminary result, suppose you allocate $$n$$ balls randomly to $$m$$ bins, and let $$K_n$$ denote the number of occupied bins. The probability of getting $$K_n=k$$ occupied bins is:

\begin{equation} \begin{aligned} p(k) &= \mathbb{P}(K_n = k) = \frac{(m)_k \cdot S(n,k)}{m^n}, \\[6pt] \end{aligned} \end{equation}

where $$(m)_k = m(m-1)(m-2) \cdots (m-k+1)$$ are the falling factorials and $$S(n,k)$$ are the Stirling numbers of the second kind. Now, to solve your problem, define the quantity:

$$T \equiv T(m_*) \equiv \min \{ n \in \mathbb{N} | K_n = m_* \}.$$

This is the number of balls that need to be randomly allocated to get $$m_*$$ occupied bins (using the stopping rule specified above). For all $$1 \leqslant m_* \leqslant m$$ we have the distribution function:

\begin{equation} \begin{aligned} F_T(t) = \mathbb{P}(T \leqslant t) &= \mathbb{P}(K_t \geqslant m_*) \\[6pt] &= \sum_{k=m_*}^{m} \frac{(m)_k \cdot S(t,k)}{m^t}. \\[6pt] \end{aligned} \end{equation}

Thus, we obtain the mass function:

\begin{equation} \begin{aligned} p_T(t) = \mathbb{P}(T = t) &= \mathbb{P}(T \leqslant t) - \mathbb{P}(T \leqslant t-1) \\[6pt] &= \sum_{k=m_*}^{m} \Bigg[ \frac{(m)_k \cdot S(t,k)}{m^t} - \frac{(m)_k \cdot S(t-1,k)}{m^{t-1}} \Bigg] \\[6pt] &= \sum_{k=m_*}^{m} \frac{(m)_k}{m^t} \Bigg[ S(t,k) - m \cdot S(t-1,k) \Bigg] \\[6pt] &= \sum_{k=m_*}^{m} \frac{(m)_k}{m^t} \Bigg[ \frac{1}{k!} \sum_{i=0}^{k} (-1)^i {k \choose i} \Big( (k-i)^{t} - m \cdot (k-i)^{t-1} \Big) \Bigg] \\[6pt] &= \sum_{k=m_*}^{m} \frac{(m)_k}{m^t} \Bigg[ \frac{1}{k!} \sum_{i=0}^{k} (-1)^i {k \choose i} (k-i)^{t-1} (k-m-i) \Bigg] \\[6pt] &= \frac{1}{m^t} \sum_{k=m_*}^{m} \frac{(m)_k}{k!} \sum_{i=0}^{k} (-1)^i {k \choose i} (k-i)^{t-1} (k-m-i) \\[6pt] &= \frac{1}{m^t} \sum_{k=m_*}^{m} \sum_{i=0}^{k} (-1)^i \frac{1}{i!} \frac{(m)_k}{(k)_i} (k-i)^{t-1} (k-m-i). \\[6pt] \end{aligned} \end{equation}

This is the probability mass function for the random variable $$T$$, which represents the number of balls allocated until you have $$1 \leqslant m_* \leqslant m$$ occupied bins.

Computing this distribution: This distribution involves the Stirling numbers of the second kind, so it involves some computational challenges. You can create a function to compute the log-probabilities for the distribution using some special R packages. The following code gives a function to compute the log-CDF and log-PDF of the distribution up to a specified upper-bound.

#Load required libraries
library(matrixStats);
library(gmp);
library(VGAM);
library(ggplot2);

#Create function for log-Stirling numbers
LS2  <- function(t, k) {
if (t < k) { -Inf } else { log(gmp::Stirling2(t,k)) } }

#Create function for log-difference
logdiff <- function(l1, l2) { l1 + VGAM::log1mexp(l1-l2); }

#Create function for distribution of T
DIST <- function(mstar, m, tupper) {

#Check that argument values are valid
if (!is.numeric(m))     { stop('Error: m must be numeric') }
if (m < 0)              { stop('Error: m must be non-negative') }
if (!is.numeric(mstar)) { stop('Error: mstar must be numeric') }
if (mstar < 0)          { stop('Error: mstar must be non-negative') }
if (mstar > m)          { stop('Error: mstar must not be larger than m') }

#Compute matrix of terms
LLL <- matrix(-Inf, nrow = tupper, ncol = m);
MMM <- rep(-Inf, tupper);
PPP <- rep(-Inf, tupper);
for (t in 1:tupper) {
for (k in mstar:m)  {
LLL[t, k] <- LS2(t,k) - t*log(m) + lgamma(m+1) - lgamma(m-k+1); }
MMM[t] <- matrixStats::logSumExp(LLL[t, ]);
if (t == 1) { PPP[t] <- MMM[t] } else {
if (MMM[t] >= MMM[t-1]) { PPP[t] <- logdiff(MMM[t], MMM[t-1]); } } }

#Give function outputs
list(log.cdf = MMM, log.pdf = PPP); }


We can implement an example for the particular values $$m=20$$ and $$m_* = 12$$ showing all argument values $$t=1,2,...,60$$. We generate a bar-plot of the probability mass function in this case. It should be possible to compare this distribution with a Monte Carlo simulation to confirm its correctness.

#Set parameter values
m      <- 20;
mstar  <- 12;
tupper <- 40;

#Generate the distribution
DDD  <- DIST(mstar, m, tupper);
DATA <- data.frame(T = 1:tupper, logprob = DDD$log.pdf); #Plot the mass function THEME <- theme(plot.title = element_text(hjust = 0.5, size = 14, face = 'bold'), plot.subtitle = element_text(hjust = 0.5, face = 'bold')); ggplot(aes(x = T, y = exp(logprob)), data = DATA) + geom_bar(stat = 'identity', fill = 'red') + THEME + ggtitle('Probability Mass Function') + xlab('T') + ylab('Probability'); • I believe the question makes sense when you consider the filling of the buckets with stopping rule that M buckets are non-zero, and then consider the number of balls needed,$T$, as the variable. Then your function for$p(m|t,n)$, when integrated (summed) to become a CDF, solves the problem. That is because the CDF for$T$is related to the CDF of$M$. $$P(M < m|t,n) = P(T > t |m,n)$$ the probability that you need more than$t$balls to fill$m$out of$n$buckets equals the probability that with$t$balls you fill less than$m$buckets. – Sextus Empiricus Sep 1 '19 at 21:26 Here is a result using simulations in R, using parameters N=40 and M=20 N=40 M=20 n=1e5 res=replicate(n,{ k=h=0 #k=total balls, h=unique balls cnt=cnt2=rep(0,N) while (h<=M) { tmp=sample(1:N,1) cnt2[tmp]=1 if (sum(cnt2[cnt2!=0])>M) { break } else { cnt[tmp]=cnt[tmp]+1 k=sum(cnt[cnt!=0]) h=sum(cnt2[cnt2!=0]) } } return(k) }) • This answer assumes that buckets are filled until the M+1 bucket would have been filled, which may be what the OP is asking. – bjschoenfeld Aug 31 '19 at 18:27 This problem is ill-posed, as is. In some sense, one could say that given M, the probability that the number of balls is K, i.e. $$p(K|M) = 0$$. Here's why. Assume N is fixed. Consider Bayes Theorem: $$p(K|M) = \frac{p(M|K)p(K)}{p(M)} = \frac{p(M|K)p(K)}{\Sigma_K p(M|K)}$$. $$p(M|K)$$ is related to the Birthday Problem: given $$K$$ people/balls, what is probability that there will be $$M$$ distinct birthdays/buckets with at least one ball? This can be computed, see here. $$p(M|K) = \frac{\binom{N}{M}{K \brace M}N!}{N^K}$$, where $$\binom{N}{M}$$ is the binomial coefficient and $${K \brace M}$$ is the stirling number of the second kind. Clearly $$K >= M$$, but there is no upperbound on $$K$$. No matter how balls $$K$$ we start with, $$M$$ could be 1, 2, ..., or N. This means that $$\Sigma_K p(M|K)$$ is an infinite sum over $$K$$ of positive values $$p(M|K)$$, with 2 cases: 1) it diverges, i.e. is infinite or 2) it converges and has a finite positive value. In case 1, we are dividing by an infinite value, so loosely speaking, $$p(K|M) = 0$$. In case 2, we need to consider the final term in Bayes Theorem: $$p(K)$$, i.e. what is the probability we started with $$K$$ balls? Given the problem, we don't know. We could say that any $$K$$ is as likely as any other. Since $$K >= M$$, there are infinitely many possible values, so $$p(K) = 0$$ (and we should probably start using measure theory). Thus $$p(K|M) = \frac{finitevalue*0}{otherfinitevalue} = 0$$. Now, if we decide on some $$p(K)$$, like $$K$$ is any positive integer less than 100 or less than $$N$$ and are all equally likely, then all the terms on the right-hand side of Bayes Theorem are finite and we could determine $$p(K|M)$$. If $$K <= 100$$ and all are equally likely, then $$p(K) = \frac{1}{100}$$ and $$p(K|M) = \frac{\frac{\binom{N}{M}{K \brace M}N!}{N^K}*\frac{1}{100}}{\Sigma_K p(M|K)} = \frac{\frac{\binom{N}{M}{K \brace M}N!}{N^K}*\frac{1}{100}}{\Sigma_K \frac{\binom{N}{M}{K \brace M}N!}{N^K}}$$. • The statement following "this means that" is mistaken. What matters is the likelihood of the data, not the sum of likelihoods over the entire hypothesis space. The conclusion that$p(K\mid M)=0$is obviously incorrect because for any$M$the distribution of$K\$ is discrete. – whuber Aug 30 '19 at 12:09