Distribution of number of non-zero counts of a multinomial distributed set The multinomial count distribution for a set of $M$ categories is
$P((n_1, \dots,n_M)|N)=\left(\frac{N!}{n_1!\dots n_m!}\prod_{i=1}^Mp_i^{n_i}\right)\delta\left( \sum_{i=1}^M n_i-N\right)\;,$
where category $i=1, \dots, m$ appears $n_i$ times, and the Dirac $\delta$-function appears to enforce the sum constraint that $\sum_{i=1}^M n_i=N$.
What is the distribution, $P(m)$ of the number of non-zero counts, $m=\sum_{i=1}^M\mathrm{sgn}[n_i]$, where $\mathrm{sgn}(x)=1$ for $x>0$ and $0$ for $x=0$.
I couldn't find anything in a quick internet search, but I am still hoping there may be a closed form solution, even for the case where the $p_i$ are all distinct.
One approach is to focus instead on the distribution, $P_0(m_0)$, of zero counts, $m_0$, from which I think $P(m)=1-P_0(M-m)$. 
A specific subset of $m_0$ categories is denoted $\alpha=\{\alpha_1, \dots, \alpha_{m_0}\}$, where $\alpha_i\in\{1, \dots, M\}$. There are $\binom{M}{m_0}$ distinct $\alpha$ subsets. For a given $\alpha$, the probability to obtain one of the $\alpha$ categories in a single count sample is $p_\alpha=\sum_{i=1}^{m_0}p_{\alpha_i}$. Thus, the probability to obtain something other than the $\alpha$ categories across $N$ samples is $(1-p_\alpha)^N$. This representation ignores permutations of the $\alpha$ subset. There are $m_0!$ such permutations. 
If someone knows a solution, please share. In lieu of a solution, I could use help in the combinatorics in the sum of $(1-p_\alpha)^N$ over subsets $\alpha$ and their orderings.
 A: I don't think this question has a closed-form solution in general. However, there is one special case where $p_i=p=\frac{1}{M}, \forall p$, e.g., a fair 6-faces die. In such case, we can use the stars and bars. Assuming $m <= min(M, N)$.
Say we have M balls to put into N buckets. No ball should be left behind. We have totally $\binom{N + M -1}{N}$ assignments. This is the denominator. Having exactly $m$ non-empty buckets can be viewed as a two-step process. First, select $m$ buckets, $\binom{M}{m}$ possible selections. And then arrange the $N$ ball(s) in the selected buckets, $\binom{N-1}{m-1}$. For example, if $M=5, N=3$, there will be following cases:

*

*Only one bucket is non-empty. You can select any of the three buckets and put all balls in. #possibilities = 3

*Only two buckets are non-empty. Selecting two buckets has 3 possibilties. Putting 5 balls in 2 buckets where no empty bucket has 4 possibilties. #possibilities = 12

*All three buckets are non-empty. Just arrange 5 balls in three buckets. Again, no empty buckets. #possibilities = 6.

The probablities are then, e.g., $p(m=3) = 6 / \binom{5+3-1}{5} = 0.2857$.
R code for this. Note since now, I am using $k$ instead $M$.
multinomial_non0_count_v1 <- function(k, n) {
  stop = min(k, n)
  total = choose(n + k -1, n)
  for (m in 1:stop) {
    prob = choose(k, m) * choose(n-1, m-1) / total
    print(paste(m, prob))
  }
}

The above function is not numerically stable for large k,n,m because the binomial coefficients become too large to compute. However, we can calculate the probability by expanding the the binomial coefficients and canceling through division. $p = \frac{[k,...,m+1][k-1,...,m][n,...,n-m+1]}{[n+k-1,...,n][k-m,...,1]}$. Abusing the notation of $[x,...,y]$, here it means $x(x-1)(x-2)...y$, the product of continuous integers between $x$ and $y$, where, $x >= y$.
prob <- function(m, k, n) {
  stopifnot(m <= min(k,n))
  if (m == k) {
    res = exp(sum(log(seq(n, n-m+1))) - sum(log(seq(n+k-1, n))))
  } else {
    minuend = sum(log(seq(k, m+1))) + sum(log(seq(k-1, m))) + sum(log(seq(n, n-m + 1)))
    subtrahend = sum(log(seq(n+k-1, n))) + sum(log(seq(k-m, 1))) 
    res = exp(minuend - subtrahend)
  }
  res
}

