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$)
 A: (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

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