Hay un método algebraico simple y elegante. Se reduce a poco más que el emparejamiento y reemplazo de patrones repetidos (con patrones muy simples), lo que lo hace eficiente (al menos para problemas pequeños como este).
El objetivo es calcular la función generadora para el número de caras que aparecen $k$ o más veces en $n$ lanzamientos de un "dado" generalizado que tiene un número finito (digamos $m$) de caras con probabilidades $\mathbf{p}=(p_1, p_2, \ldots, p_m)$ de aparecer. Esto, por definición, es un polinomio
$$g(u; n,k,\mathbf{p}) = g_0 + g_1u + \cdots + g_ju^j + \cdots$$
donde $g_j$ es la probabilidad de exactamente $j$ caras apareciendo $k$ o más veces.
La pregunta particular concierne el caso $k=3$, $n=10$, y $\mathbf{p}=\left(\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6}\right)$. Para ilustrar el proceso, consideremos un "dado" de dos caras con "caras" indexadas por $1$ y $2$ teniendo probabilidades $p=\left(\frac{1}{2},\frac{1}{2}\right)$ de aparecer: una "moneda justa." Disminuyamos $n$ a $5$ y $k$ a $2$ para mantener las expresiones cortas.
Comencemos con una función generadora del dado en sí,
$$f(\mathbf{x}; \mathbf{p}) = p_1 x_1 + p_2 x_2 + \cdots +p_m x_m.$$
Esto es un polinomio en las $m$ variables $\mathbf{x}=\left(x_1, \ldots, x_m\right)$. Para la moneda justa,
$$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.$$
Podemos leer todos los resultados posibles de $n$ lanzamientos expandiendo $f(\mathbf{x},\mathbf{p})$ a la potencia $n$-ésima.
$$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}.$$
Esto indica que hay una probabilidad de $\frac{1}{32}$ de ver cinco unos, una probabilidad de $\frac{5}{32}$ de ver cuatro unos y un dos, y así sucesivamente.
El truco es evitar hacer este cálculo completo, que es puro forcejeo. Notemos que no necesitamos la mayoría: solo necesitamos reconocer los términos donde la potencia de $x_1$ o $x_2$ es $k=2$ o mayor. El resto lo podemos desechar. Una forma de lograr esto--que se hace riguroso por el concepto algebraico de un ideal polinómico--es introducir $m$ variables nuevas, digamos $u_1, u_2, \ldots, u_m$. Simplemente declaramos que $u_i = x_i^k$ para cada $i$ y usamos esto para simplificar la expansión. Además, las reglas de álgebra implican que podemos realizar la simplificación calculando $f(\mathbf{x},\mathbf{p})^n$ de forma recursiva como
$$f(\mathbf{x},\mathbf{p})^n = f(\mathbf{x},\mathbf{p})f(\mathbf{x},\mathbf{p})^{n-1}$$
y haciendo la simplificación en cada paso. En el ejemplo, este procedimiento da como resultado
$$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.$$
Podemos pulir esto con algunas conversiones más. Primero, la aparición de cualquier potencia de $u_i$ igual a $2$ o mayor simplemente refleja el hecho de que $x_i^k$ se encontró varias veces. Por ejemplo, $u_1^2$ significa que $x_1^k$ apareció dos veces--es decir, $x_1$ apareció al menos $2k$ veces. Además, después de haber calculado $f\left(\mathbf{x},\mathbf{p}\right)^n$ podemos ignorar cualquier potencia de los $x_i$ que aún esté presente. Esto se hace fácilmente estableciendo todos los $x_i=1$. En el ejemplo, el resultado de estas dos operaciones es
$$\frac{5 u_2 u_1}{8}+\frac{3 u_1}{16}+\frac{3 u_2}{16}.$$
Todo lo que queda por hacer es seguir la cantidad total de cada monomio: $u_2u_1$ tiene grado $2$, lo que indica que dos caras separadas aparecieron al menos $k$ veces. $u_1$ y $u_2$ tienen grado $1$, lo que indica que exactamente una cara apareció al menos $k$ veces. El grado es fácil de calcular: establece todas las variables $u_i=u$ para una única variable $u$. En el ejemplo, hemos encontrado así que
$$g(u;n,k,\mathbf{p}) = \frac{5 u^2}{8}+\frac{3 u}{8}.$$
(Interpretación: después de cinco lanzamientos de una moneda justa, hay una probabilidad de $3/8$ de que exactamente una cara se observe dos o más veces y una probabilidad de $5/8$ de que ambas caras se observen dos o más veces. Estos resultados se pueden verificar fácilmente.)
Para el caso de la pregunta, con diez lanzamientos de un dado justo, estos cálculos dan como resultado
$$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}.$$
A partir de esto, podemos calcular fácilmente cualquier probabilidad que deseemos. Por ejemplo, el evento en el que dos caras distintas aparecen cada una tres o más veces (y todas las demás caras dos veces o menos) tiene una probabilidad igual al coeficiente de $u^2$ y, por lo tanto, es $52415/139968 \approx 0.374478.$
Usando un sistema algebraico computacional, el código es notablemente corto. La siguiente es una solución en Mathematica. Todo el trabajo se realiza en las tres líneas centrales, siguiendo exactamente la descripción anterior.
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)
La ejecución tomó 0.2 segundos.