En primer lugar, permite calcular $f(m)$, el número de maneras de manejar los dos naipes a cada uno de $m$ de las personas que son de todos los pares. Podemos hacer esto usando un equipo de la siguiente manera. Si además vamos a $f(m,r)$ ser la probabilidad de que todos los $m$ de las personas tienen las parejas cuando las cartas de una baraja de $4r$ tarjetas, a continuación,
$$
f(m,r) = r\binom42\big(f(m-1,r-1)+(m-1)f(m-2,r-1)\big)
$$
El primer sumando representa caminos, donde el primer jugador coge un par que nadie más tiene, y el segundo cuentas por caminos por los que algún otro jugador también recibe ese par. La anterior ecuación recursiva permite calcular $f(m)=f(m,13)$ rápidamente a través de la programación dinámica.
A continuación, utilice la generalizada del principio de inclusión exclusión. Dejando $E_i$ ser el caso de que el $i^{th}$ jugador tiene un par, entonces
\begin{align}
P(n,k)
&=\sum_{i=0}^n(-1)^{i-k}\binom{i}{k}\binom{n}{i}\mathbb P(E_1\cap E_2\cap \dots \cap E_i)\\
&=\sum_{i=0}^n(-1)^{i-k}\binom{i}{k}\binom{n}{i}\frac{f(i)}{\frac{52!}{2^i(52-2i)!}}
\end{align}
Usted puede ver esto en acción aquí, simplemente haga clic en el botón ejecutar en la parte superior e introduzca el valor deseado de $n$.
También, existe la siguiente "forma cerrada" para $P(n,k)$:
$$
\boxed{P(n,k) = \binom{n}{k}\sum_{i=k}^n(-1)^{i-k}\binom{n-k}{i-k}\frac{2^i(52-2i)!}{52!}\cdot i!\cdot [x^i](1+6x+3x^2)^{13}}
$$
Aquí, la notación $[x^i]f(x)$ significa que el coeficiente de $x^i$ en el polinomio $f(x)$.