Picking pairs of socks from a drawer. There are $n$ socks in a drawer, of $m$ different colours. Initially, the probability of picking a sock of colour $c_i$ at random is $\mathbb{P}(c_i) \cdot 2r$ socks are picked at random, without replacement. 
What is the probability, $\mathbb{P}(pairs)$, that $r$ pairs of socks are picked? (i.e. all socks are paired)
If $r = 1$ (2 socks are chosen):
$\mathbb{P}(pairs)$ = $\mathbb{P}$(2 are $c_1$ or 2 are $c_2$ or ... or 2 are $c_m$) ~ $\sum_{k}^{m} {\left[\mathbb{P}(c_k)\right]^2}$
How can this be generalized for r?

In short, although I can write a solution to your problem, I think that simplifying it may be quite intractible. 
Say you choose $i_j$ socks of colour $j$. Then $\sum_{j=1}^Mi_j=2r$, $i_j\leq n_j$ for all $j$ and the probability that you get pairs in this particular way is $\left(\begin{array}{c} 2r \\ i_1\;i_2\;\ldots\;i_m\end{array}\right)\frac{n_1(n_1-1)\ldots(n_1-i_1+1)}{n^{i_1}}\frac{n_2(n_2-1)\ldots(n_2-i_2+1)}{n^{i_2}}\ldots\frac{n_m(n_m-1)\ldots(n_m-i_m+1)}{n^{i_m}}$, where in the case $i_j=0$ the corresponding term is equal to 1. The multinomial coefficient is there because you can choose the socks in different orders. Therefore $\mathbb{P}[\text{Pairs}]=\sum_{(i_j)_j:\begin{array}{c} \sum_{j=1}^Mi_j=2r \\ i_j\leq n_j \\ i_j \text{ even}\end{array}}\left(\begin{array}{c} 2r \\ i_1\;i_2\;\ldots\;i_m\end{array}\right)\frac{n_1(n_1-1)\ldots(n_1-i_1+1)n_2(n_2-1)\ldots(n_2-i_2+1)}{n^{i_2}}\ldots\frac{n_m(n_m-1)\ldots(n_m-i_m+1)}{n^{i_m}}=\\\sum_{(i_j)_j:\begin{array}{c} \sum_{j=1}^Mi_j=2r \\ i_j\leq n_j \\ i_j \text{ even}\end{array}}\left(\begin{array}{c} 2r \\ i_1\;i_2\;\ldots\;i_m\end{array}\right)\frac{n_1(n_1-1)\ldots(n_1-i_1+1)n_2(n_2-1)\ldots(n_2-i_2+1)\ldots n_m(n_m-1)\ldots(n_m-i_m+1)}{n^{2r}}$
To be clear the summation notation means we are summing over all vectors $(i_1,\ldots,i_M)$ that satisfy the given constraints. In the limit as $n_j\rightarrow\infty$ for all $j$ but $r$ stays constant we get the approximation $\sum_{(i_j)_j:\sum_{j=1}^Mi_j=2r, i_j \text{ even}}\left(\begin{array}{c} 2r \\ i_1\;i_2\;\ldots\;i_m\end{array}\right)\frac{n_1^{i_1}n_2^{i_2}\ldots n_m^{i_m}}{n^{2r}}$. 
This is well studied, it is the probability that a multinomial distribution has all terms even. As far as I know there is no clean answer to this, which means I doubt there is a clear simplification of the exact answer, although it can be calculated in the case of specific values of $r$, $n_i$. 
