Suppose I roll 10 dice (six-sided). What is the probability of seeing k faces that appear more than (or equal to) 3 times?
For example, if I see 1, 1, 1, 2, 2, 2, 3, 4, 5, 6 then I have observed 2 faces that appeared more than (or equal to) 3 times (1 and 2).
I guess I can try running a simulation but was wondering if there is a neat way to calculate this analytically.
There is a simple elegant algebraic method. It amounts to little more than repeated pattern matching and replacement (with very simple patterns), making it efficient (at least for small problems like this one).
The aim is to compute the generating function for the number of faces appearing $k$ or more times out of $n$ rolls of a generalized "die" that has a finite number (say $m$) of faces with probabilities $\mathbf{p}=(p_1, p_2, \ldots, p_m)$ of appearing. This, by definition, is a polynomial
$$g(u; n,k,\mathbf{p}) = g_0 + g_1u + \cdots + g_ju^j + \cdots$$
where $g_j$ is the chance of exactly $j$ faces appearing $k$ or more times.
The particular question concerns the case $k=3$, $n=10$, and $\mathbf{p}=\left(\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6}\right)$. To illustrate the process, though, consider a two-sided "die" with "faces" indexed by $1$ and $2$ having probabilities $p=\left(\frac{1}{2},\frac{1}{2}\right)$ of appearing: a "fair coin." Let's decrease $n$ to $5$ and $k$ to $2$ to keep the expressions short.
Begin with a generating function of the die itself,
$$f(\mathbf{x}; \mathbf{p}) = p_1 x_1 + p_2 x_2 + \cdots +p_m x_m.$$
This is a polynomial in the $m$ variables $\mathbf{x}=\left(x_1, \ldots, x_m\right)$. For the fair coin,
$$f\left(\left(x_1,x_2\right); \left(\frac{1}{2},\frac{1}{2}\right)\right) = \frac{1}{2}x_1 + \frac{1}{2}x_2.$$
We can read all possible outcomes of $n$ throws by expanding $f(\mathbf{x},\mathbf{p})$ to the $n^\text{th}$ power.
$$f(\mathbf{x},\mathbf{p})^5 = \frac{x_1^5}{32}+\frac{5}{32} x_2 x_1^4+\frac{5}{16} x_2^2 x_1^3+\frac{5}{16} x_2^3 x_1^2+\frac{5}{32} x_2^4 x_1+\frac{x_2^5}{32}.$$
This says there is a $\frac{1}{32}$ chance of seeing five ones, a $\frac{5}{32}$ chance of four ones and one two, and so on.
The trick is to avoid doing this full calculation, which is pure brute force. Notice that we don't need most of it: we only need to recognize the terms where the power of $x_1$ or $x_2$ is $k=2$ or greater. The rest we can throw away. One way to accomplish this--which is made rigorous by the algebraic concept of a polynomial ideal--is to introduce $m$ new variables, say $u_1, u_2, \ldots, u_m$. We simply declare that $u_i = x_i^k$ for each $i$ and use this to simplify the expansion. Moreover, the rules of algebra imply we may perform the simplification by computing $f(\mathbf{x},\mathbf{p})^n$ recursively as
$$f(\mathbf{x},\mathbf{p})^n = f(\mathbf{x},\mathbf{p})f(\mathbf{x},\mathbf{p})^{n-1}$$
and doing the simplification at each step. In the example this procedure yields
$$f(\mathbf{x},\mathbf{p})^5 = \frac{1}{32} u_1^2 x_1+\frac{5}{32} u_1^2 x_2+\frac{5}{16} u_2 u_1 x_1+\frac{5}{16} u_2 u_1 x_2+\frac{5}{32} u_2^2 x_1+\frac{1}{32} u_2^2 x_2.$$
We may polish this off with a few more conversions. First, the appearance of any power of $u_i$ of $2$ or greater just reflects the fact that $x_i^k$ was found multiple times. For instance, $u_1^2$ means $x_1^k$ appeared twice--that is, $x_1$ appeared at least $2k$ times. Furthermore, after having computed $f\left(\mathbf{x},\mathbf{p}\right)^n$ we can ignore any powers of the $x_i$ still lying around. That is easily done by setting all the $x_i=1$. In the example the result of these two operations is
$$\frac{5 u_2 u_1}{8}+\frac{3 u_1}{16}+\frac{3 u_2}{16}.$$
All that remains to do is track the total degree of each monomial: $u_2u_1$ has degree $2$, indicating that two separate faces appeared at least $k$ times. $u_1$ and $u_2$ have degree $1$, indicating that exactly one face appeared at least $k$ times. The degree is easy to compute: set all the variables $u_i=u$ for a single variable $u$. In the example we have thereby found that
$$g(u;n,k,\mathbf{p}) = \frac{5 u^2}{8}+\frac{3 u}{8}.$$
(Interpretation: after five flips of a fair coin, there is a $3/8$ chance that exactly one face will be observed two or more times and a $5/8$ chance that both faces will be observed two or more times. These results are easily checked.)
For the setting in the question, with ten tosses of a fair die, these calculations yield
$$g\left(u; 10, 3, \left(\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6}\right)\right)=\frac{4375 u^3}{209952}+\frac{52415 u^2}{139968}+\frac{225559 u}{419904}+\frac{175}{2592}.$$
From this we may easily compute any probabilities we wish. For instance, the event where two distinct faces each appear three or more times (and all other faces two times or fewer) has probability equal to the coefficient of $u^2$ and therefore is $52415/139968 \approx 0.374478.$
Using a computer algebra system, the code is remarkably short. The following is a Mathematica solution. All the work is done in the three middle lines, exactly following the foregoing description.
g[y_, p_List, n_Integer, k_Integer] /; n >= 1 && k >= 1 :=
Module[{f, x, u, a, i},
f = Sum[p[[i]] Subscript[x, i], {i, 1, Length[p]}]/Total[p];
Nest[(Expand[f #] /. {Power[Subscript[x,i_], a_] /; a == k -> Subscript[u,i]}) &, 1, n]
/. {Subscript[x, _] -> 1, Power[Subscript[u, _], _] -> y, Subscript[u, _] -> y}
]
g[u, ConstantArray[1/6, 6], 10, 3]
Execution took 0.2 seconds.
