Solution (source code)

= Solution

The totals are <independent> binomials, $n\sim\operatorname{Bin}(2N,q)$ and $m\sim\operatorname{Bin}(2M,p)$. Thus
$$
L(p,q)=\binom{2N}n\binom{2M}m q^n(1-q)^{2N-n}p^m(1-p)^{2M-m}.
$$
Differentiating the log <likelihood>, or maximizing the two factors separately, gives
$$
\boxed{\hat p=\frac m{2M},\qquad\hat q=\frac n{2N}.}
$$
This remains the constrained <maximum> at zero or one when a total is at an endpoint. At interior parameters the sample information, in coordinate order $(p,q)$, is
$$
\boxed{I(p,q)=\begin{pmatrix}2M/[p(1-p)]&0\\0&2N/[q(1-q)]\end{pmatrix}.}
$$
For example $-\mathbb E\ell_{pp}=\mathbb E[m/p^2+(2M-m)/(1-p)^2]=2M/[p(1-p)]$; the cross <derivative> is zero. As both sample sizes grow with nondegenerate interior probabilities, the joint approximation is normal with mean $(p,q)$ and <covariance> $\operatorname{diag}(p(1-p)/(2M),q(1-q)/(2N))$. The <independence> of the two estimators is exact, inherited from the two data samples.