How to calculate probability of $N$ pairings between two lists of random integers? Say I generate two lists $L_1, L_2$ of length $l_1, l_2$ of integers uniformly selected from $[1,m]$.
I then count the pairings of same values between the lists by crossing out an entry in $L_1$ that matches an entry in $L_2$ (and also crossing out that entry in $L_2$) so that those particular entries are eliminated from further consideration.
I'd like to calculate, given $l_1, l_2,$ and $m$ the distribution of the number of pairings.
Simulating this is trivial, but I have no idea how to approach this mathematically (and apologize for having no work to show in that regard).
Is there a way to do this directly?
 A: We shall represent the lists as multisets using indicator vectors. For instance, the lists $L_1=[1,1,2]$ and $L_2=[1,2,2]$, for $m = 5$, will be represented as $$
L_1=(2,1,0,0,0)\\
L_2=(1,2,0,0,0) \enspace.
$$
Take the entry-wise minium $(1,1,0,0,0)$
which sums to 2, the number of pairings for this instance. 
Using this representation we can model the random variable $L_1$ (resp. $L_2$) as a multinomial random variable with $l_1$ trials (resp. $l_2$) and uniform probability vector $\boldsymbol p=(1,1,\ldots,1)/m$ of length $m$, with probability mass function $f_1(\boldsymbol x) = \binom{l_1}{\boldsymbol x} /m^{l_1}$ (resp. $f_2(\boldsymbol x) = \binom{l_2}{\boldsymbol x}/ m^{l_2}$), in which $\binom{l_1}{\boldsymbol x}$ for a vector $\boldsymbol x$ of length $m$ is the multinomial coefficient $l_1!/\boldsymbol x! = l_1!/x_1!x_2!\cdots x_m!$. 
The probability of observing $N = i$ pairings, for $i\in\{0,1,\ldots,\min(l_1,l_2)\}$, is thus:
$$ \sum_{\substack{\min(\boldsymbol x,\boldsymbol y) = \boldsymbol z }} f_1(\boldsymbol x) f_2(\boldsymbol y)  \enspace,
$$
in which the sum is over $\boldsymbol x\in C(l_1,m),\boldsymbol y\in C(l_2,m),$ and $\boldsymbol z\in C(i,m)$, in which $C(i;m)$ is the set of all $m$-weak compositions of $i$. 
The expression can be expanded into
$$
\frac{l_1! l_2!}{m^{l_1+l_2}}
\sum_{\substack{\min(\boldsymbol x,\boldsymbol y) = \boldsymbol z }}
\frac{1}{\boldsymbol x!\boldsymbol y!} \enspace.
$$
In this form it may be computationally intensive to calculate, but perhaps it may be possible to simplify the expression further into something that can be evaluated efficiently. Here's Mathematica code to compute it:
Needs["Combinatorica`"]
m^-(l1 + l2) l1! l2!
Total[Function[z, 
   Total[(1/Times @@ (Join @@ #1!) &) /@ 
     Select[Tuples[{Compositions[l1, m], Compositions[l2, m]}], 
      z == (Min @@ #1 &) /@ Thread[#1] &]]] /@ Compositions[i, m]]

