Two contests, an extension of the Coupon Collector's problem Coupon Collector's Problem
Let $X$ be the number of coupons drawn with replacement from an urn containing $N$ distinct coupons until each coupon has been drawn at least once, winning the coupon contest at the store.
It's clear to me that the expectation of $X$ is $N H_N$ where $H_N$ is the $N$-th harmonic number.
Consider now a slightly different problem:
Let $X$ and $Y$ be the numbers of coupons drawn with replacement from two different urns containing $N$ and $M$ distinct coupons (respectively) until every coupon has been drawn at least once from their respective urns, winning both coupon contests at the store. Coupons are drawn two at a time, one from each contest.

The expectation of $\max(X, Y)$ is what I'm looking for, the expectation of how many double-draws were made to win both contests. The difficulty I'm having is that both contests are differently distributed, the "max formula" they taught in class those many semesters ago doesn't fly here since they're not identically distributed. Even harder still, is the fact that, though I understand the result of the coupon collector's problem, I'm not familiar enough in its construction to confidently tackle this problem.
Anyone have some immediate hints or even a full solution to this? I'm still trying struggling to get $P(X = x)$ and $P(Y = y)$.
 A: We start by describing the distributions of $X$ and $Y$ separately.  Write $X=\sum_{i=1}^N X_i,$ where $X_i$ is the time to collect the $i$-th coupon.  The $X_i$'s are independent geometrically distributed: that is, $X_i\sim G(p_i)$ where $p_i=(N+1-i)/N$. We have $\mathbf{E}X_i=1/p_i$ and $Var X_i=(1-p_i)/p_i^2$, so
$$\mathbf{E}X = \sum_{i=1}^N \frac{1}{p_i} = \sum_{i=1}^N \frac Ni = N H_1(N)$$
and
$$Var X = \sum_{i=1}^N \frac{1-p_i}{p_i^2} = \sum_{i=1}^N \left(\frac{N^2}{i^2} - \frac Ni\right) = N^2 H_2(N)-N H_1(N)$$
where $H_k(N)$ denotes the "order-$k$" harmonic number.  Of course the analogous formula holds for $Y$.
The reason this problem is complex is that the pdf's of $X$ and $Y$ are convolutions, the cdf's are sums of these, and the cdf for $\max(X,Y)$ depends on the cdf's of $X$ and $Y$:
$$\mathrm{Pr}(\max(X,Y)\le k)=\mathrm{Pr}(X\le k)\,\mathrm{Pr}(Y\le k).$$
This is unlikely to have a closed form.  However, for large $M$ and $N$ you can verify that the sequences $\{X_i\}$ and $\{Y_j\}$ satisfy the conditions of the Lyapunov central limit theorem, so $X$ and $Y$ can be approximated by normal distributions with the mean and variance indicated above.  This may suffice depending on how detailed an analysis you need.
If you really want to pursue an exact expression, the following approach may work.  The Laplace transform $L_X(z)=\mathbf{E}z^X$ captures complete information about the distribution of $X$; in particular $\mathrm{Pr}(X=k)$ (resp. $\mathrm{Pr}(X\le k)$) is the coefficient of $z^k$ in $L_X(z)$ (resp. $L_X(z)/(1-z)$).  In this case, we have $\mathbf{E}z^{X_i}=p_i z/(1-q_i z)$ where $q_i=1-p_i$, so
$$L_X(z)=\prod_{i=1}^N \mathbf{E}z^{X_i}
= \frac{z^N\prod p_i}{\prod (1-q_i z)}.$$
However, you're going to need to analyze the Hadamard (i.e. termwise) product of $L_X(z)/(1-z)$ and $L_Y(z)/(1-z)$, rather than the ordinary product, so any analysis is likely to require complex integration (see, e.g. Formal power series coefficient multiplication.)
