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

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

-
Just too note that this probability is 1 whenever $n\geq mk$. and you may have problems with values of $n$ which are close to $mk$. certain approximations may break down. –  probabilityislogic Feb 13 '12 at 3:37
I could be wrong, but isnt the answer you want basically the cdf of a multinomial distribution? if so you are probably looking at the "multivariate" incomplete bet function –  probabilityislogic Feb 13 '12 at 3:51
@probabilityislogic: We can make a slightly stronger statement: $n>m(k-1) \Rightarrow p(n,m,k)=1$; also $n<m \Rightarrow p(n,m,k)=0$. Yes the problem involves a multinomial distribution, but I fail to see how "at least one of the values is drawn at least $k$ times" translates into something computable as a CDF. –  fgrieu Feb 13 '12 at 11:22
At least one is the "reciprical" of none drawn k times or more. Thus we have our probability as one minus a multivariate cdf. See @daniel's answer. –  probabilityislogic Feb 13 '12 at 15:39

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$.

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.

-
I admire the attempt to produce an efficient exact solution. The crux of this is your recurrence: where does this come from? It does not appear to obtain correct answers. For instance, my simulation with 10^7 iterations (reported in another comment) has a standard error of 0.00015 and your result is (0.64732 - 0.64597) / 0.00015 = 8.9 standard errors too high. Thus, either the simulation is badly wrong or your solution is incorrect. (I suspect you may be double-counting.) What, for instance, do you get for the values of $W(4,4,3,k)$? They should be 54, 78, 81 for $k=2,3,4$ respectively. –  whuber Feb 13 '12 at 15:31
@whuber I am more than willing to admit I am wrong, but where do you get $W(4,4,3,2)=54$? This is the number of ways of placing 4 balls in 3 bins such that no bin has two or more balls. I cannot find even one way. –  deinst Feb 13 '12 at 16:08
@whuber I believe you mean $k=3,4,5$. In which case, yes I get those numbers. –  deinst Feb 13 '12 at 16:11
Sorry, my interpretation of $k$ was off by 1. Since you get those numbers, it looks like your approach is basically right, at least for these small values. But to tell whether it is, we really need a clear explanation of the recurrence. –  whuber Feb 13 '12 at 16:58
(+1) In light of your improved explanation, I have to agree. But before that explanation was available, there wasn't any good basis to decide which to favor. I'll take a well-reasoned theoretical argument any time over simulation results :-). –  whuber Feb 13 '12 at 20:23
show 1 more comment

(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

-
(Updated to reflect corrected code.) A simulation with nIterations=10^7 yields 0.647 23 (that one took 40 minutes over lunch :-). The first 3.5 digits of this result should be correct, strongly indicating the Poisson approximation (0.645 9) is actually more accurate than the Binomial (0.639 9). –  whuber Feb 8 '12 at 17:13
Excel gives 0.640.., which appears significantly different from 0.645.. given by the simulation. Perhaps the hypothesis made for the Excel formula is not so correct after all. –  fgrieu Feb 8 '12 at 17:44
Please don't trust Excel for such extended calculations. It is known to have problems with statistical distribution functions, especially for extreme values of their arguments. (I was shocked to find that it returned anywhere near a reasonable value...) R (1 - pbinom(19, 3000, 1/300)^300) and Mathematica (1 - CDF[BinomialDistribution[3000, 1/300], 19]^300 // N) both give 0.639877. –  whuber Feb 8 '12 at 19:17
I do not trust Excel too much either, but here it gives 0.639877005.. which agrees to 10 figures with Mathematica 0.63987700518.. thus the problem is either with the approximation made or (less likely) the simulation. –  fgrieu Feb 8 '12 at 20:14
Yes, it's interesting that Excel is getting a good answer. Clearly the Binomial approximation is decent but not entirely accurate. At least in this case, the Poisson approximation is highly accurate. It's good to mistrust theoretical calculations and simulations, even when they agree, so I invite you to produce a simulation of your own as a check! (You could do this with an Excel macro, but that's a bit of a pain...) –  whuber Feb 8 '12 at 20:33

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$$

-
I believe the exponent should be an m, not a k as you have. I have included another step that makes this more clear. I may be missing something, or misunderstanding your question, so please clarify if you still don't agree. –  Daniel Johnson Feb 8 '12 at 16:40
To help resolve these differences, I have edited the question to show the value Excel actually returns for the formula as given (with $k$ in the exponent). Putting $m$ in the exponent yields a value of 0.6398; compare this to the simulation results. –  whuber Feb 8 '12 at 16:56
So the rest of your question is how to represent this expression in R? I don't use R, so I can't help you with that, but i would guess there is a binomial cdf function that you can use and if not, the summation formulation should be pretty trivial to code. –  Daniel Johnson Feb 8 '12 at 19:50
@DanielJohnson: Thanks for your explanation and correction. I can now use R for simulation (as well as computing the theoretical approximation), thanks to whuber. The only part remaining is how not to depend on an approximation that is not so good (it introduce a relative error on p of nearly 1%, and likely much worse in other setups). –  fgrieu Feb 8 '12 at 20:21

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.

-
I have added code that performs the normal approximation from the paper at the end of my answer (it was too big to fit here). I get $p(3000, 300, 20)\approx 0.6468276$ (an error of about 0.0005). –  deinst Feb 14 '12 at 17:35