Consider a 'balls-and-bins model' where we allocate $n$ balls at random to $m$ bins and each ball has a probability $\theta$ of occupying its allocated bin (and a corresponding probability of $1-\theta$ of 'falling through' its bin so that it does not occupy that bin). Note that this is the model used in the 'extended occupancy problem'.
Question: What is the distribution of the maximum number of balls occupying any single bin, and how do we compute its probabilities?
To formalise this model, let $U_1,...,U_n \sim \text{IID U} \{ 1,...,m \}$ denote the random allocations of the $n$ balls to $m$ bins, and let $Q_1,...,Q_n \sim \text{IID Bern}(\theta)$ be the corresponding indicators that a ball occupies its allocated bin (i.e., it does not fall through the bin). The counts for the number of balls occupying each bin (and the number of balls that fall through bins) are given by:
$$\begin{align} N_1 &= \sum_{i=1}^n Q_1 \cdot \mathbb{I}(U_i = 1), \\[6pt] &\quad \quad \quad \quad \quad \vdots \\[6pt] N_m &= \sum_{i=1}^n Q_m \cdot \mathbb{I}(U_i = m), \\[6pt] N_* &= n - \sum_{i=1}^n Q_i. \\[6pt] \end{align}$$
The vector of counts (including the count $N_*$ for the balls falling through bins) follows the distribution:
$$\mathbf{N} = (N_1,...,N_m, N_*) \sim \text{Mu} \Bigg( n, \Big( \frac{\theta}{m}, ..., \frac{\theta}{m}, 1-\theta \Big) \Bigg),$$
and the maximum count of balls occupying a bin is $M_n \equiv \max (N_1,...,N_m)$. To facilitate our nalysis, we will refer to the distribution of this statistic as the MaxCount distribution and we will denote its probability mass function as follows:
$$\mathbb{P}(M_n = s) = \text{MaxCount}(s|n, m, \theta).$$
In the sections below we will show a mixture formulation for this distribution and then show how to compute probabilities from this distribution. This relates closely to the distribution of the maximum-count of a uniform-multinomial distribution, which is examined in Bonetti, Cirillo and Ogay (2019). We will establish computation for this distribution by using a mixture representation in terms of this latter distribution, and using the iterative method in that paper to compute the probabilities.
A mixture formulation for the MaxCount distribution: In the 'classical case' where $\theta = 1$ all of the balls occupy their bins with probability one, and so we have $N_*=0$ almost surely (i.e., no balls fall through their bins). In this case the statistic $M_n$ corresponds to the maximum-count of a uniform-multinomial distribution. In the broader case, the more general distribution can be written as a mixture of the distribution for the classical case.
$$\text{MaxCount}(s|n, m, \theta) = \sum_{n_\text{eff} = s}^n \text{MaxCount}(s|n_\text{eff}, m, 1) \cdot \text{Bin}(n_\text{eff} | n, 1-\theta).$$
The intuition behind this mixture representation is simple --- the value $n_\text{eff}$ denotes the effective number of balls that occupy their allocated bins, and so the probability for the maximum-count is obtained by applying the law of total probability conditional on the effective number of balls. The mixture representation allows us to compute the general MaxCount distribution from the classical MaxCount distribution.
Computing the MaxCount distribution: The classical case for the MaxCount distribution (i.e., when $\theta = 1$) corresponds to the distribution of the maximum-count of a uniform-multinomial distribution. An iterative method for computing the cumulative distribution function for this distribution is given in Bonetti, Cirillo and Ogay (2019) (pp. 6-7). We can use this iterative method to compute the CDF for the classical case (which fortunately computes an array of probabilities for all values of $n$ up to our specified value) and then we can apply the above mixture representation to take account of the parameter $\theta$ in the distribution.
The MaxCount distribution is available in the occupancy package in R using standard syntax for probability functions (dmaxcount, pmaxcount, qmaxcount and rmaxcount). The distribution uses the parameters size for $n$, space for $m$, and prob for $\theta$. The algorithm in the function does all computation in log-space, and computes the log-probabilities by combining the iterative algorithm in Bonetti, Cirillo and Ogay (2019) with the mixture representation above.
Here is an example where we compute the probabities for an urn model with $n=20$ balls and $m=6$ bins with $\theta = 1$. As can be seen from the plot, it is highly likely that the maximum-count in the bins will be between 4-8 balls in this case.
#Set the parameters
n <- 20
m <- 6
theta <- 1
#Compute the probabilities from the MaxCount distribution
PROBS <- dmaxcount(0:n, size = n, space = m, prob = theta)
barplot(PROBS, names.arg = 0:n, col = 'blue',
xlab = 'Maximum Count', ylab = 'Probability')
