No creo que exista una solución general de forma cerrada. Se puede encontrar la solución exacta mediante programación dinámica. Usted puede aproximar la distribución mediante el cálculo de la media y la varianza y el uso de una aproximación normal. Y, por último, para el caso en que $p_1 = p_2 = \ldots p_M = 1/M$ se puede encontrar una solución de forma cerrada utilizando números de Stirling del segundo tipo.
Solución exacta mediante programación dinámica
Defina $S_k = \sum_{i=1}^k n_i$ , $I_k = \sum_{i=1}^k \mathbb{I}(n_i > 0)$ . Mediante programación dinámica calcularemos iterativamente la probabilidad de que el $k$ las primeras variables suman $n$ y que $m$ de la primera $k$ son distintos de cero,
$$ d(k, n, m) = P(S_k = n, I_k = m) $$
Inicializar la tabla de programación dinámica con,
$$ d(0, n, m) = \begin{cases} 1 & n=0, m=0 \\ 0 & \text{otherwise} \end{cases} $$
A continuación, recurse sobre $k=1,\ldots,M$ ,
$$ d(k, n, m) = \operatorname{binomial}(0; N-n, \tilde{p}_k)d(k-1, n, m) + \sum_{r=1}^{n}\operatorname{binomial}(r; N-n+r,\tilde{p}_k)d(k-1, n-r, m-1), \quad m \leq k $$
donde $\tilde{p}_{k} = p_{k} / \sum_{i=k+1}^M p_i$ . El primer término tiene en cuenta el caso en el que la variable k'th es cero y, por tanto, no contribuye con nuevas variables distintas de cero, y el segundo término (la suma) tiene en cuenta todos los casos en los que la variable k'th es distinta de cero. A continuación se muestra una implementación ingenua en R. Los bucles se pueden reorganizar para que el bucle final termine en $r$ en lugar de $m$ para evitar computar las mismas densidades binomiales múltiples veces, aún así no creo que esto se pueda hacer mejor que $\mathcal{O}(M^2N^2)$ .
p <- c(0.1,0.1,0.2,0.2,0.4)
N <- 6
M <- 5
p_tilde <- p / rev(cumsum(rev(p)))
d <- array(0, dim = c(M+1,N+1,M+1))
d[1,1,1] = 1
for (k in 1:M) {
for (n in 0:N) {
for (m in 0:k) {
d[k+1, n+1, m+1] <- dbinom(0, N-n, prob = p_tilde[k]) * d[k, n+1, m+1]
if (m > 0 && n > 0) {
r <- 1:n
d[k+1, n+1, m+1] <- d[k+1, n+1, m+1] +
sum(dbinom(r, r+N-n, prob = p_tilde[k]) * d[k, n-r+1, m])
}
}
}
}
Verificar la solución mediante simulación
> table(colSums(rmultinom(1e6, size = N, prob = p) > 0 )) / 1e6
1 2 3 4 5
0.004245 0.114157 0.449484 0.374784 0.057330
> d[M+1,N+1,]
[1] 0.000000 0.004226 0.114734 0.449280 0.374160 0.057600
Solución de aproximación
Para utilizar una aproximación normal debemos calcular la media y la varianza. En primer lugar, la media se puede calcular como,
$$ \mathbb{E}[m]= \mathbb{E}\left[\sum_{i=1}^M\mathbb{I}(n_i>0)\right] = \sum_{i=1}^M\mathbb{E}[\mathbb{I}(n_i>0)] = \sum_{i=1}^M 1-(1-p_i)^N. $$
A continuación, la varianza, \begin{equation} \begin{split} \mathbb{V}[m] &= \mathbb{V}\left[\sum_{i=1}^M\mathbb{I}(n_i>0)\right] \\ &= \sum_{i=1}^M\mathbb{E}[\mathbb{I}(n_i>0)] \\ &= \sum_{i=1}^M\mathbb{V}[\mathbb{I}(n_i>0)] + 2 \sum_{1\leq i < j \leq M} \operatorname{Cov}(\mathbb{I}(n_i>0), \mathbb{I}(n_j>0)) \\ &= \sum_{i=1}^M (1-(1-p_i)^N)(1-p_i)^N + 2 \sum_{1\leq i < j \leq M} \operatorname{Cov}(1-\mathbb{I}(n_i>0), 1-\mathbb{I}(n_j>0)) \\ &= \sum_{i=1}^M (1-(1-p_i)^N)(1-p_i)^N + 2 \sum_{1\leq i < j \leq M} \left[(1-p_i-p_j)^N - (1-p_i)^N(1-p_j)^N\right] \end{split} \end{equation}
Se trata de un $\mathcal{O}(N^2)$ algoritmo. Si $M$ es muy grande, puede omitir los términos de covarianza, que en general debería esperar que sean pequeños, y puede calcular la aproximación en $\mathcal{O}(N)$ .
computed_variance <- sum((1-p)^N*(1-(1-p)^N))
for (i in 1:(M-1)) {
for (j in (i+1):M) {
computed_variance <- computed_variance + 2*((1-p[i]-p[j])^N-(1-p[i])^N*(1-p[j])^N)
}
}
> computed_variance
[1] 0.602115
> var(colSums(rmultinom(1e6, size = N, prob = p) > 0 ))
[1] 0.602461
Estuche para uniformes
Por último, en el caso uniforme $p_i=p\,\forall i$ debemos contar de cuántas maneras podemos distribuir $N$ objetos numerados en $M$ número utilizando exactamente $m$ de los contenedores. Esto se puede hacer en un procedimiento de dos pasos, primero elija $m$ de la $M$ bins, esto se puede hacer en $M \choose m$ formas. A continuación, distribuya el $N$ en cada uno de los $m$ con al menos un objeto en cada recipiente, esto se puede hacer en $m!S(N,m)$ donde $S(N,m)$ es el número de Stirling del segundo tipo. Obtenemos la siguiente fórmula
$$ p(m) = {M \choose m} m! S(N,m) / M^N $$
Ejemplo, $N=4$ , $M=5$ , $m=3$ ( $S(4,3) = 6$ ):
> table(colSums(rmultinom(1e6, size = 4, prob = rep(1/5, 5)) > 0)) / 1e6
1 2 3 4
0.008013 0.224771 0.575393 0.191823
> choose(5, 3) * factorial(3) * 6 / 5^4
[1] 0.576
