When sampling a categorical distribution with $k$ outcomes of equal likelihood $n$ times, what is the probability of sampling all categories?

Question

Suppose I have a multinomial distribution with $k$ categories and $n$ trials where the probabilities of each category are the same ($1/k$). I would like to know the probability that every category has at least one trial, i.e. the sum of multinomial probabilities where $x_i\geq 1$ for all $i=1,...,k$. For $n<k$ this probability is zero and the case where $n=k$ this is straightforwardly $n!/n^n$, but I am having trouble generalizing this to greater $n$ values.

Though it likely follows from the answer to the previous question, I am also curious about finding the likelihood that all categories are sampled at least $m$ times, i.e. $x_i\geq m$ for all $i=1,...,k$.

Answer 1

We fix $k$ and let $p_{n,i}$ denote the probability that after $n$ draws, we have drawn exactly $i$ unique categories. Thus, what we are interested in is $p_{n,k}$. We will calculate $p_{n,i}$ dynamically, filling a matrix with $n$ rows and $k$ columns.

For the first row, we have $p_{1,1}=1$ and $p_{1,i}=0$ for $2\leq i\leq k$.

For rows below the first ($n\geq 2$), we first consider $i=1$. We can only have drawn $1$ category after $n$ draws if we have had exactly $1$ category after $n-1$ draws and drawn the same category again, with a probability of $\frac{1}{k}$. Thus, $$ p_{n,1} = p_{n-1,1}\times\frac{1}{k}. $$

For $i\geq 2$, we can arrive at having drawn exactly $i$ categories after $n$ draws in one of two mutually exclusve ways:

1. We have already drawn $i$ categories after $n-1$ draws, and re-drew one of these in the $n$-th draw, for a total probability of $p_{n-1,i}\times\frac{i}{k}$.
2. We have drawn $i-1$ categories after $n-1$ draws, and drew one of the previously undrawn categories in the $n$-th draw, for a total probability of $p_{n-1,i-1}\times\frac{k-(i-1)}{k}$.

Overall, $$ p_{n,i} = p_{n-1,i}\times\frac{i}{k} + p_{n-1,i-1}\times\frac{k-(i-1)}{k}. $$

This is easily calculated, and a simulation gives results that match closely. In R:

kk <- 5
nn <- 10

probs <- matrix(0,nrow=nn,ncol=kk,dimnames=list(1:nn,1:kk))
probs[1,1] <- 1

for ( draw in 2:nn ) {
    probs[draw,1] <- probs[draw-1,1]*1/kk
    for ( n_unique in 2:kk ) {
        probs[draw,n_unique] <- probs[draw-1,n_unique]*n_unique/kk +
            probs[draw-1,n_unique-1]*(kk-(n_unique-1))/kk
    }
}

set.seed(1)
sims <- replicate(1e5,length(unique(sample(1:kk,nn,replace=TRUE))))

probs[nn,]
table(factor(sims,levels=1:kk))/length(sims)

Result:

> probs[nn,]
          1           2           3           4           5 
0.000000512 0.001046528 0.057323520 0.419082240 0.522547200 
> table(factor(sims,levels=1:kk))/length(sims)

      1       2       3       4       5 
0.00000 0.00102 0.05739 0.42157 0.52002

Essentially, what we are doing is a Markov Chain with $n$ steps through a state space with $k$ possible states, the $i$-th state meaning "we have drawn $i$ unique categories". Something similar would probably be possible for your more general question about having seen each category $m$ times, but there, we would need to keep track of how often we have already drawn each category, so we have to look at far more possible states.

Answer 2

Occupancy problem

This is a special case of the occupancy problem.

Start with $k$ empty urns, and place $n$ times a ball into a randomly chosen urn, what is the probability that $x$ urns are non-empty?

You case is the special case $x=k$, the probability that all are non-empty.

---

Computation with Stirling numbers of second kind

This allocation of balls can be considered by

• First partitioning the $n$ balls into $k$ unlabelled groups.

  The number of ways this can be done is the Stirling number of the second kind.

• Second assign $k$ labels to the $k$ unlabelled groups.

  The number of ways this can be done is the factorial $k!$.

For example, you draw $10$ balls, let's number them $1, 2, \dots , 9, 10$, and place them into three unlabeled categories. One possible partition would be $\{1,4,5\} - \{2,6,7,9,10\} -\{3,8\}$. This could correspond to the particular way that the balls get distributed into the $3$ classes. The Stirling number gives you all the possible ways to do this. You multiply by $3! = 6$ because there are $6$ ways to assign the labels $1,2,3$ to the $3$ unlabeled groups in the partition.

The total number of ways to have non-zero classes is

$$n_{\text{non-zero}} = k! {{n}\brace{k}} $$

---

Solution, the expression for the probability

Then your probability becomes

$$P(\text{non-zero}) = \frac{n_{\text{non-zero}}}{n_{\text{total}}} = \frac{k!{{n}\brace{k}}}{k^n}$$

the computation can be done by this summation

$$P(\text{non-zero}) = \sum_{i=0}^k {{k}\choose{i}} (-1)^i \left( \frac{k-i}{k} \right)^n$$

Approximation

We could use an approximation like this

$$k! {{k+t}\brace{k}} \sim \frac{k^{k+2t} \sqrt{2 \pi k}}{e^k \cdot 2^t \cdot t!} \left(1 + \frac{2t^2}{3k} \right)$$

and

$$P(\text{non-zero}) \sim \frac{(k/2)^{n-k} \sqrt{2 \pi k}}{e^k \cdot (n-k)!} \left(1 + \frac{2(n-k)^2}{3k} \right) $$

Alternative approximation

We could also keep on adding balls into the categories until all $k$ categories are filled, and consider the distribution of how many balls $n$ it takes to fill all $k$ categories. See also the coupon collector's problem.

The probability that it takes $n$ or fewer steps is equal to the probability that after $n$ steps all categories are full.

The waiting time to advance the number of non-zero categories by one step, follows the geometric distribution with parameter $p = x/k$ where $x$ is the number of categories that are still open. The average waiting time is $1/p$ and the variance is $(1-p)/p^2$.

The total waiting time is sum of all those geometric distributed waiting times and can be approximated with a Gaussian distribution. The mean and variance are:

$$\begin{array} {} \mu &=& \sum_{x=1}^k (x/k)^{-1} &=& k H_k\\ \sigma^2 &=& \sum_{x=1}^k (1-x/k)(x/k)^{-2} &=& k^2 H_k^{(2)}- kH_k \end{array}$$

where $H_k$ is the $k-th$ harmonic number and $H_k^{(2)}$ is the generalized harmonic number.

Then we approximate it with a Gaussian distribution and look for the probability that the variable with this mean and variance is below $n$. So we look for the CDF of standardized deviation from the mean of $\frac{n - \mu }{\sigma}$ which has probability $$P(\text{non-zero}) \approx \Phi\left( \frac{n - \mu}{\sigma} \right) = \Phi\left( \frac{n/k - H_k}{\sqrt{H_k^{(2)}- H_k/k}} \right)$$ and we can approximate those harmonic numbers by $H_k^{(2)} \approx \pi^2/6$ and $H_k \approx log(k) + \gamma + \frac{1}{2k}$ then

$$P(\text{non-zero}) \approx \Phi\left( \frac{n/k - (\gamma + \log(k) + 0.5/k)}{\sqrt{\pi^2/6 - (\gamma + \log(k) + 0.5/k)/k }} \right)$$

Edit: improved alternative approximation

It has been shown by several that the tail of the coupon collector's problem approaches a Gumbel distribution. So we could try to make the same approach as above but use a Gumbel distribution where we match the mean and variance by

$$\begin{array} {} \mu &=& k H_k &\approx & k (\gamma + ln(k) + 0.5/k)\\ \sigma^2 &=& k^2 H_k^{(2)}- kH_k & \approx & k^2 \frac{\pi^2}{6}- \mu \end{array}$$

giving (with $\beta = \sqrt{6\sigma^2/\pi^2}$)

$$P(\text{non-zero}) = \frac{k!{{n}\brace{k}}}{k^n} \approx e^{-e^{-(x-\mu+\beta\gamma)/\beta}}$$

The simulation below shows that this works very well.

This also gives an interesting (new?) approach in approximating Stirling numbers of the second kind.

$${{n}\brace{k}} \approx \frac{k^n}{k!} e^{-e^{-(x-\mu+\beta\gamma)/\beta}}$$

Simulation

The simulation below shows that the direct formula works (as expected).

The same computation using an estimate of the Stirling numbers of the second kind does not work very well (I wonder if I made a mistake somewhere, or maybe there just needs to be some more terms added to the estimate).

To my surprise, the estimate with the normal distribution is not that bad (especially when compared to the estimate with the estimated Stirling numbers of the second kind)

### function to compute stirling number
Stirling2 = function(n,k) {
  i = 0:k
  terms <- (-1)^i * choose(k,i) * (k-i)^n
  sum(terms)/factorial(k)
} 
Stirling2 = Vectorize(Stirling2)


### function to compute stirling number with terms from prob included
### this reduces the size of the numbers and allows to compute larger numbers
### multiplied with factorial(k)
### devided by k^n
Stirling3 = function(n,k) {
  i = 0:k
  terms <- (-1)^i * choose(k,i) * ((k-i)/k)^n
  sum(terms)
} 
Stirling3 = Vectorize(Stirling3)



### function to compute harmonic
harmonic = function(k,a) {
  sum(1/(1:k)^a)
}

### function to simulate a draw
Sim <- function(n,k) {
  x = sample(1:k, n, replace = TRUE)
  cats = c(1:k) %in% x
  prod(cats)
}

### function to compute probability 
SimP <- function(n,k,m) {
  sum(replicate(m,Sim(n,k)))/m
}
SimP = Vectorize(SimP, vectorize.args = c("n"))


### settings
set.seed(1)
k = 40       
n = k+(0:100)*3

### simulutions
psim = SimP(n,k,10^3)

### compute with Stirling
pcomp = Stirling2(n,k)*factorial(k)/k^n
pcomp = Stirling3(n,k)


### estimate with Stirling
logpest = (n-k)*log(k/2)+0.5*log(2*pi*k) - k - lfactorial(n-k)
pest = exp(logpest) * (1+2*(n-k)^2/3/k)

### estimate with normal
mu = k*harmonic(k,1)
var = k^2*harmonic(k,2) - mu
pest2 = pnorm((n-mu)/var^0.5)

std_x = (n/k - (0.57721 + log(k) + 0.5/k))/
        sqrt(pi^2/6 - (0.57721 + log(k) + 0.5/k)/k)
pest3 = pnorm(std_x)

### estimate with gumbel distribution
mu = k*(0.57721 + log(k) + 0.5/k)
var = k^2*pi^2/6 - mu
beta = sqrt(var*6/pi^2)
loc = mu - beta * 0.57721
pest4 = VGAM::pgumbel(n, location = loc , scale = beta)

### plot
plot(n,psim, log = "")
lines(n,pcomp)
lines(n,pest, lty = 2)
lines(n,pest3, lty = 3)
lines(n,pest4, lty = 1, col = 3)
# lines(n,pest2, lty = 3, col = 2) # plots the same as pest3
title("k = 40 and n from 40 to 340")

legend(40,1, c("simulated values", "computation with stirling numbers", "estimate with normal distribution", "estimate with stirling numbers approximation"),
       cex = 0.85, pch = c(1,NA,NA,NA), lty = c(NA,1,3,2))