M buckets out of N have at least one ball

Question

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?

Answer 1

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');

Answer 2

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