$n$ players are each dealt two cards — what's the probability that $k$ of them have a pair? Each of $n\leq26$ players is dealt $2$ cards from a standard $52$-card poker deck. What is $\textrm{P}\left(n,k\right)$, the probability that exactly $k$ of the $n$ players have a pair?
(A pair is a hand like $8 \clubsuit, 8 \heartsuit$ or $K \clubsuit, K \diamondsuit$.)

This question was previously asked at Poker.SE without a satisfactory answer: https://poker.stackexchange.com/questions/4087/probability-of-x-pocket-pairs-at-a-table-of-n-people-nlhe
 A: The probability that a particular player is dealt a pair is $\frac{4-1}{52-1}=\frac{1}{17}$ so the expected number of pairs is $\frac{n}{17}$ 
I suspect the distribution is not going to be far away from Binomial with parameters $n$ and $p=\frac{1}{17}$. This is not quite correct (for example when $n=26$ the probability of exactly $25$ pairs is $0$ rather than $4.24\times 10^{−30}$) but I would be surprised if it were substantially misleading
As an example, here are a couple of simulations of $100000$ million deals to $6$ players 
Pairs                    0        1         2         3         4         5         6 
Simulation 1            69459    26065     4094      367       13         2         0
Simulation 2            69675    25918     4063      329       15         0         0
Binomial probability  0.695066 0.260650 0.040727 0.003394 0.000159 0.0000040 0.00000004

and, given the uncertainty of simulation, these do not look noticeably different.  The five pocket pairs out of six seen at the link in your link was indeed unlikely
A: First, let us compute $f(m)$, the number of ways to deal two cards to each of $m$ people which are all pairs. We can do this using a computer as follows. If we further let $f(m,r)$ be the probability that all $m$ people have pairs when dealt from a deck of $4r$ cards, then
$$
f(m,r) = r\binom42\big(f(m-1,r-1)+(m-1)f(m-2,r-1)\big)
$$
The first summand accounts for ways where the first player gets a pair that no one else has, and the second accounts for ways where some other player also gets that pair. The above recursive equation allows you to compute $f(m)=f(m,13)$ quickly using dynamic programming. 
Next, use the generalized principle of inclusion exclusion. Letting $E_i$ be the event that the $i^{th}$ player has a pair, then
\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}
You can see this in action here, just click the run button at the top and enter your desired value of $n$.

Also, there is the following "closed form" for $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}}
$$
Here, the notation $[x^i]f(x)$ means the coefficient of $x^i$ in the polynomial $f(x)$. 
