How do you find quantiles in this balls-in-bins problem?

Question

I need to calculate the expected number of hash collisions with a range for a software project. I think this is a reformulation of the birthday problem, as follows.

Suppose you have $n$ balls allocated at random to $d$ bins. What is the greatest value $m$ such that the probability of getting at least $m$ balls in the same bin exceeds one-half?

Answer 1

In order to be as responsive as possible to the applied scenario of interest to you, I am going to frame my analysis in terms of data objects that take on a random hash value, so that objects with the same hash value give hash collisions. The analysis is no different if you use balls/bins instead of data objects/hash values, so please interpret with whatever applied context you wish.

Your problem here involves an analysis of the maximum count from a uniform-multinomial distribution; it can be framed mathematically as follows. Let $X_1,...,X_n \sim \text{IID U} \{ 1,...,d \}$ be a set of $n$ uniform draws from $d$ hash values. Define the count values $N_r = \sum_{i=1}^n \mathbb{I}(X_i = r)$ for each hash value and observe that the count vector has a uniform-multinomial distribution:

$$\mathbf{N} = (N_1,...,N_d) \sim \text{Mu} \Big( n, (\tfrac{1}{d},...,\tfrac{1}{d}) \Big).$$

Define the largest count over a single hash value as $M_n \equiv \max (N_1,...,N_d)$. You are looking for the largest value $m$ such that $\mathbb{P}(M_n \geqslant m) \geqslant 0.5$, which is a quantile of the distribution of $M_n$.

Unfortunately, the distribution of the uniform-multinomial maximum-count is a complicated distribution that does not have a simple closed form. An iterative method for computing the distribution of the maximum is examined in Bonetti, Cirillo and Ogay (2019) (pp. 6-7). It is programmed in R in the occupancy package using standard syntax for probability functions (see related question here). Here I will give an example with $n=20$ and $d = 6$, but you can alter these numbers if you prefer.

#Set the parameters
n <- 20
d <- 6

#Compute and plot the cumulative probabilities from the MaxCount distribution
library(occupancy)
CUMPROBS <- pmaxcount(0:n, size = n, space = d)
names(CUMPROBS ) <- sprintf('M[%s]', 0:n)
plot(0:n, CUMPROBS, col = 'blue', pch = 20, cex = 1.5,
     xlab = 'Maximum Count', ylab = 'Cumulative Probability')
abline(a = 0.5, b = 0, lty = 2, lwd = 2)

#Compute quantile
QUANT <- qmaxcount(0.5, size = n, space = d)
QUANT

[1] 6

In order to get $\mathbb{P}(M_n \geqslant m) \geqslant 0.5$ we can solve the equivalent problem $\mathbb{P}(M_n \leqslant m-1) \leqslant 0.5$. As you can see from the above plot of the cumulative distribution function, the largest value of $m$ that satisfies this requirement is $m=6$. In other words, if we have $n=20$ pieces of data distributed randomly over $d=6$ distinct hash values, there is at least a 50% chance that the maximum number of pieces of data with the same hash is at least six.

Answer 2

If you are willing to work with estimations of the probabilities there is a rather simple solution with no need of calculation. Here is some R code you can play with.

birthdays = function(n, d){ # Asigns bds to n people
  sample(1:d, n, replace = TRUE )
}

collisions = function(bds, m){ #Check if there are m or more people who share birthda
  any(table(bds)>= m)
}

estimate_prob =function(n,d,m) { 
  mean( replicate(2000, collisions(birthdays(n,d), m)) ) # we generate 2000 samples of bds anc calculate proportion of those where ther are m or more people sharing bday. You can change number o samples to a greater number for more acuracy in estimation.
}

The code is in R. The function birthdays samples a birthday in a year of d days for a total of n people. Collisions returns TRUE if and only if m or more share the same bday. The function estimate_pro estimates the probability of m colissions between n people and a d days year. You can check that for $d = 365$ and $n = 23$ the probability of sharing two people sharing birthday is aproximately $0.5$.

estimate_prob(23, 365, 2)

Then you can use this function and for given $n$ and $d$ estimate the probability of m collisions on a grid. Here is an example

prob = rep(0,10)

for( m in 1:10) {
  prob[m] = estimate_prob(100, 365, m)
}
plot(1:10, prob, xlab = "m")
abline(0.5, 0, col = "red")

This gives the following plot which shows that for $d = 365$ and $n = 100$ the probability of three or more collisions is about $0.6$ and four or more is around $0.05$