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 Aug 30, 2019 at 2:35
• If N=M, then I think there is no solution. Might as well be infinite number of balls. Aug 30, 2019 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... Aug 30, 2019 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. Aug 30, 2019 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, 2019 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{aligned} p(k) &= \mathbb{P}(K_n = k) = \frac{(m)_k \cdot S(n,k)}{m^n}, \\[6pt] \end{aligned}

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{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}

Thus, we obtain the mass function:

\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}

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. Sep 1, 2019 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. Aug 31, 2019 at 18:27