When sampling a categorical distribution with $k$ outcomes of equal likelihood $n$ times, what is the probability of sampling all categories? 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$.

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

