This question is more tricky than it sounds. To solve it, I will use here a combination of both manual methods and automated methods, in particular computer algebra tools [the mathStatica
package (of which I am an author) for Mathematica and the latter itself].
If I may change the notation slightly:
The Problem
Let $X_1$ ~ $\text{Binomial}(n,p)$ and $X_2$ ~ $\text{Binomial}(m,q)$ be independent.
Find the pmf of $|X_1-X_2|$
Given: Due to independence, the joint pmf of $(X_1, X_2)$, say $f(x_1,x_2)$, is:
Solution
Let $Y=X_1-X_2$ and $Z=X_2$. Then, the joint pmf of $(Y,Z)$, say $g(y,z)$, is:
where Transform
is a mathStatica
function that derives the joint pmf using the Method of Transformations. Deriving the domain of support of $Y$ and $Z$ is a bit more tricky. To make things clearer, here is a rough plot that illustrates the (smoothed continuous version of) the domain of support:
This suggests two cases:
The density of $Y=X_1-X_2$ is then obtained by summing out $Z$ in each part of the domain:
Finally, to find the pmf of $|Y|$, the pmf for strictly positive values will be:
and when $Y=0$:
Summary
The pmf of $|X_1-X_2|$, say $\phi(y)$ is:
with domain of support $Y$ = {0, 1, ..., max$(m,n)$}.
All done.
Monte Carlo check
It is always a good idea to check ones work using Monte Carlo methods. Here, for instance, are 100,000 pseudo-random drawings from each of $X_1$ and $X_2$, given some parameter assumptions:
x1data = RandomVariate[BinomialDistribution[12, .1], 100000];
x2data = RandomVariate[BinomialDistribution[ 7, .9], 100000];
Next, compare the empirical distribution of $|X_1-X_2|$ (red triangles) to the theoretical density $\phi(y)$ (blue dots) derived above, given the same parameter assumptions:
Looks good :)