Odds of drawing at least k identical values among m after n draws?

Question

We draw $n$ values, each equiprobably among $m$ distinct values. What are the odds $p(n,m,k)$ that at least one of the values is drawn at least $k$ times? e.g. for $n=3000$, $m=300$, $k=20$.

Note: I was passed a variant of this by a friend asking for "a statistical package usable for similar problems".

My attempt: The number of times a particular value is reached follows a binomial law with $n$ events, probability $1/m$. This is enough to get odds $q$ that a particular value is reached at least $k$ times [Excel gives $q\approx 0.00340$ with =1-BINOMDIST(20-1,3000,1/300,TRUE)]. Given that $n\gg k$, we can ignore the fact that odds of a value being reached depends on the outcome for other values, and get an approximation of $p$ as $1-(1-q)^m$ [Excel gives $p\approx 0.640$ with =1-BINOMDIST(20-1,3000,1/300,TRUE)^300].

update: the exponent was wrong in the above, that's now fixed

Is this correct? (now solved, yes, but the approximation made leads to an error in the order of 1% with the example parameters)

What methods can work for arbitrary parameters $(n,m,k)$? Is this function available in R or other package, or how could we construct it? (now solved, both exactly for moderate parameters, and theoretically for huge parameters)

I see how to do a simulation in C, what would be an example of a similar simulation in R? (now solved, a corrected simulation in R and another in Python gives $p\approx 0.647$)

Answer 1

There are almost certainly easier ways, but one way of computing the value precisely is compute the number of ways of placing $n$ labeled balls in $m$ labeled bins such that no bin contains $k$ or more balls. We can compute this using a simple recurrence. Let $W(n,j,m',k)$ be the number of ways of placing exactly $j$ of the $n$ labeled balls in $m'$ of the $m$ labeled bins. Then the number we seek is $W(n,n,m,k)$. We have the following recurrence: $$W(n,j,m',k)=\sum_{i=0}^{k-1}\binom{n-j+i}{i}W(n,j-i,m'-1,k)$$ where $W(n,j,m',k)=0$ when $j<0$ and $W(n,0,0,k)=1$ as there is one way to pace no balls in no bins. This follows from the fact that there are $\binom{n-j+i}{i}$ ways to choose $i$ out of $n-j+i$ balls to put in the $m'$th bin, and there are $W(n,j-i,m'-1,k)$ ways to put $j-i$ balls in $m'-1$ bins.

The essence of this recurrence if that we can compute the number of ways of placing $j$ out of $n$ balls in $m'$ bins by looking at the number of balls placed in the $m'$th bin. If we placed $i$ balls in the $m'$th bin, then there were $j-i$ balls in the previous $m'-1$ bins, and we have already calculated the number of ways of doing that as $W(n,j-i,m'-1,k)$, and we have $\binom{n-j+i}{i}$ ways of choosing the $i$ balls to put in the $m'$th bin (there were $n-j+i$ balls left after we put $j-i$ balls in the first $m'-1$ bins, and we choose $i$ of them.) So $W(n,j,m',k)$ is just the sum over $i$ from $0$ to $k-1$ of $\binom{n-j+i}{i}W(n,j-i,m'-1,k)$.

Once we have computed $W(n,n,m,k)$ the probability that at least one bin has at least $k$ balls is $1-\frac{W(n,n,m,k)}{m^n}$.

Coding in Python because it has multiple precision arithmetic we have

import sympy
# to get the decimal approximation

#compute the binomial coefficient
def binomial(n, k):
    if k > n or k < 0:
        return 0
    if k > n / 2:
        k = n - k
    if k == 0:
        return 1
    bin = n - (k - 1)
    for i in range(k - 2, -1, -1):
        bin = bin * (n - i) / (k - i)
    return bin

#compute the number of ways that balls can be put in cells such that no
# cell contains fullbin (or more) balls.
def numways(cells, balls, fullbin):
    x = [1 if i==0 else 0 for i in range(balls + 1)]
    for j in range(cells):
        x = [sum(binomial(balls - (i - k), k) * x[i - k] if i - k >= 0 else 0
                 for k in range(fullbin))
             for i in range(balls + 1)]
    return x[balls]

x = sympy.Integer(numways(300, 3000, 20))/sympy.Integer(300**3000)
print sympy.N(1 - x, 50)

(sympy is just used to get the decimal approximation).

I get the following answer to 50 decimal places

0.64731643604975767318804860342485318214921593659347

This method would not be feasible for much larger values of $m$ and $n$.

ADDED

As there appears to be some skepticism as to the accuracy of this answer, I ran my own Monte-Carlo approximation (in C using the GSL, I used something other than R to avoid any problems that R may have provided, and avoided python because the heat death of the universe is happening any time now). In $10^7$ runs I got 6471264 hits. This seems to agree with my count, and is considerably at odds with whubers. The code for the Monte-carlo is attached.

I have finished a run of 10^8 trials and have gotten 64733136 successes for a probability of 0.64733136. I am fairly certain that things are working correctly.

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

const gsl_rng_type * T;
gsl_rng * r;

int
testrand(int cells, int balls, int limit, int runs) {
  int run;
  int count = 0;
  int *array = malloc(cells * sizeof(int));
  for (run =0; run < runs; run++) {
    int i;
    int hit = 0;
    for (i = 0; i < cells; i++) array[i] = 0;
    for (i = 0; i < balls; i++) {
      array[gsl_rng_uniform_int(r, cells)]++;
    }
    for (i = 0; i < cells; i++) {
      if (array[i] >= limit) {
    hit = 1;
    break;
      }
    }
    count += hit;
  }
  free(array);
  return count;
}

int
main (void)
{
  int i, n = 10;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (i = 0; i < n; i++) 
    {
      printf("%d\n", testrand(300, 3000, 20, 10000000));
    }

  gsl_rng_free (r);

  return 0;
}

EVEN MORE

Note: this should be a comment to probabilityislogic's answer, but it won't fit.

Reifying probabilityislogic's answer (mainly out of curiosity), this time in R because a foolish inconsistency is the hobgoblin of great minds, or something like that. This is the normal approximation from the Levin paper (the Edgeworth expansion should be straightforward, but it is more typing than I'm willing to expend)

# an implementation of the Bruce Levin article here limit is the upper limit on 
#   bin size that does not count
approxNorm <- function(balls, cells, limit) {
  # using N=s
  sp <- balls / cells
  mu <- sp * (1 - dpois(limit, sp) / ppois(limit, sp))
  sig2 <- mu - (limit - mu) * (sp - mu)
  x <- (balls - cells * mu) / sqrt(cells * sig2)
  p2 <- exp(-x^2 / 2)/sqrt(2 * pi * cells * sig2)
  p1 <- exp(ppois(limit, sp, log.p=TRUE) * cells)
  sqrt(2 * pi * balls) * p1 * p2
}

and 1 - approxNorm(3000, 300, 19) gives us $p(3000, 300, 20) \approx 0.6468276$ which is not too bad at all.

Answer 2

(Responding to the simulation question)

In R:

n <- 3000; m <- 300; k <- 20 # Problem parameters
nIterations <- 10000         # Number of iterations in the simulation
set.seed(17)                 # Make the output reproducible
#
# All the work is done in the following line.
#
t <- table(replicate(nIterations, any(tabulate(floor(runif(n, min=1, max=m+1))) >= k)))
t[["TRUE"]]/nIterations      # Convert the count to a proportion

Expected output:

[1] 0.6453

To see more deeply into what's going on, look at the detailed distribution in several experiments by executing this command several times:

table(tabulate(floor(runif(n, min=1, max=m+1))))

Typical output is

 2  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 21 
 2  5 11 23 23 32 48 37 33 21 19 21  9  6  5  1  2  1 

In this experiment, two values in the range 0..299 were observed 19 times (among 3000 independent draws) and one value was observed 21 times. You will find that most of the time, at least one value occurs 20 or more times. Because $1/m$ is small and $n$ is large, you should be seeing a Poisson distribution here. Indeed,

1 - ppois(k-1, n/m)^m

returns

[1] 0.6458719

Answer 3

This is a harder question if you don't have the $n\gg k$ and assuming that this makes them 'close enough' to independent to not affect the answer non-trivially. Lets proceed with these assumptions. Let $X_j \sim Binomial(n,\frac{1}{m})$ $\forall j = 1,..,m$.

$$P(\max_j X_j \geq k) = 1 - P(\max_j X_j < k)$$ $$ = 1 - P(X_1 < k,...,X_m < k)$$ and, assuming independence of the $m$ random variables, $$ = 1 - \prod^m_{j=1}P(X_j < k)$$ $$ = 1 - [P(X_1 < k)]^m$$ $$ = 1 - [\sum^{k-1}_{i=0} {n \choose i}(\frac{1}{m})^i(1-\frac{1}{m})^{n-i}]^m$$

or, if you have the binomial cdf function in the language you are using:

$$ = 1 - [Binomial\_cdf(k-1;n,\frac{1}{m})]^m$$

Answer 4

The collection of numbers has a multinomial distribution with $m$ categories and $n$ sample size. Letting $N_i$ be the number of times the $i$th category is chosen/repeated, we have $$(N_1,\dots,N_m)\sim multinomial\left(n;\frac{1}{m},\frac{1}{m},\dots,\frac{1}{m}\right)$$ Now leveraging off of @danieljohnson's answer the probability we are after is

$$p(n,m,k)=1-Pr(N_1<k,\dots,N_m<k)$$

i.e. if all numbers are repeated less than $k$ times, then none are repeated at least $k$ times. And "not none" is the same as "at least one" so we can take the probability away from one. This could be computed via a "brute force" approach, as the pmf we have is particularly simple:

$$p(n,m,k)=1-m^{-n}\sum_{N_1<k,\dots,N_m<k|N_1+\dots+N_m=n}{n\choose N_1\dots N_m}$$ $$=1-\frac{n!}{m^{n}}\sum_{N_1=0}^{k-1}\sum_{N_2=0}^{k-1}\dots\sum_{N_{m-1}=0}^{k-1}\frac{1}{N_1!N_2!\dots N_{m-1}!(n-N_1-N_2-\dots-N_{m-1})!}$$

The last formula is correct provided we interpret a negative factorial as $\pm\infty$ (consistent with the gamma function) which eliminates these from the summation.

On doing a quick google search came up with Bruce Levin's article. This gives a representation of the multinomial distribution as a collection of poisson random variables, with their sum being fixed. (note this might explain why @whuber has found that poisson approximation works better than binomial). Now, using the representation given in theorem 1 of the paper, we have:

$$p(n,m,k)=1-\frac{n!}{s^n\exp(-s)}\left[\prod_{j=1}^{m}Pr(X_j\leq k-1)\right]Pr(W=n)$$

Where $X_j\sim Poisson\left(\frac{s}{m}\right)$ and are independent, and $W=\sum_{j=1}^{m}Y_j$ is a sum of independent truncated poisson distributions - basically $Y_j$ is $X_j$ conditioned to be less than or equal to $k-1$. Note that we can simplify the general formula by noting that the terms in the product do not depend on the index $j$, and so is just a single poisson cdf raised to the power of $m$. Thus we have:

$$p(n,m,k)=1-\frac{n!}{s^n}\left[e_{k-1}\left(\frac{s}{m}\right)\right]^mPr(W=n)$$

Where $e_k(x)=\sum_{j=0}^{k}\frac{x^j}{j!}$ denotes the exponential sum function. note that because we have factorial an powers of potentially large numbers, numerically it will probably be better to work in terms of the logarithm of the second term, and then exponentiate back at the end of the calculation. Alternatively, we can choose the recommended $s=N$ as our algorithm parameter, and then make use of the stirling approximation to $n!$ - this is recommended in the paper and corresponds to "mean matching" of each poisson distribution with the multinomial cell (i.e. $E(X_i)=E(N_i)$). Then we get $\frac{n!}{n^n}\approx\sqrt{2\pi n}$.

The paper provides two approximations for $Pr(W=n)$ on based on normal approximation, and another based on edgeworth expansion. details are in the paper (see equation 4). Note though that his method allows for different probability parameters, so terms like $\frac{1}{t}\sum_{1}^t\sigma_l^2$ can be replaced with $\sigma_1^2$ and so on, which avoid unecessary computation. Note that we also have the mallows bounds provided in the paper - which can be used to check the accuracy of the approximations.