Probability Generating Function for The Difference of Two Binomially Distributed Random Variables?

Question

Suppose I have 2 random variables:

$X\sim \textrm{Bin}(m,p_1)$ and $Y\sim \textrm{Bin}(n,p_2).$

I want to find the distribution of $S=X-Y$ using the probability generating function ($PGF$) treating $S$ as a function of random variables. $$G_{S}(z) = G_X(z) G_Y(1/z)$$ Using the PGF for a single binomially distributed random variable:
$$G_{S}(z) = {[1-p_1+p_1 z]}^m \space {\left[1-p_2+\frac{p_2}{z}\right]}^n$$
Then to find the probability mass function PMF of the resulting distribution I use: $$p(k) = \Pr(S=k) = \frac{G_{S}^{(k)} (0)}{k!}$$ where the numerator is the $k$-th derivative of the $PGF$.

My question is how is it possible to evaluate this expression when $G_{S}^{(k)}$ contains $1/z$ terms that would be undefined at $0$ (this problem doesn't exist for sum of random variables).

(I'm aware there are other ways to find this distribution, but I just wanted to understand how to find it using the PGF).

Answer 1

Support of $ S:=X-Y$ is $[-n, m].$

So, the pgf would be

\begin{align}\mathsf P_s(t) &=\sum_{i=-n}^m p_it^i\\&= \frac{p_{-n}}{t^{n}}+\cdots+\frac{p_{-1}}{t}+p_0+p_1t+\cdots+p_mt^m\tag1\label 1\end{align}

In order to generate the probability sequence from the pgf, you can adhere to the general result:

Theorem: If a random variable $X$ takes a finite number of real values, say $\langle \alpha_i\rangle_{i=-n}^m,$ and the sequence is monotonically increasing, then $$p_{-n}=\lim_{t\to 0^+}(t^{-\alpha_{-n}}\mathsf P_s(t) ).$$

This is easy to see: From $\mathsf P_s(t) =\sum_{i=-n}^m p_it^{\alpha_i},$ we get $$t^{-\alpha_{-n}}\mathsf P_s(t) =p_{-n}+\sum_{i=-n+1}^m p_it^{\alpha_i-\alpha_{-n}}.\tag 2\label 2$$ Then take the limit.

To get $p_{-n+1}, $ now we need to differentiate $\eqref 2$ w.r.t. $t~\lfloor \alpha_{-n+1}-\alpha_{-n}\rfloor$ times and dividing both sides by $t^{\{\alpha_{-n+1}-\alpha_{-n}\}}$ and check the limiting value as $t\to 0^+.$

Rest of the coefficients follow similarly.

Coming to the present problem, from $\eqref 1,$

$$p_n=\lim_{t\to 0^+} (t^n\mathsf P_S(t)) . $$

From that very relation, we can see

$$ t^n\mathsf P_S(t) = p_{-n} +p_{-n+1}t^{-n+1+n} +\cdots.$$ Differentiate this and see what happens as $t$ tends to zero.

You can then investigate the others.

---

Reference:

$\rm [I]$ M. L. Esquível, Probability generating functions for discrete real-valued random variables, Teor. Veroyatnost. i Primenen., $2007,$ Volume $52,$ Issue 1, $129–149$ DOI: https://doi.org/10.4213/tvp8

Answer 2

You can see the probability generating function, when applied to the binomial distribution, as a way to do the checks and balances in the different ways that the heads and tails can be distributed (heads and tails, as in flipping a coin $n$ times and the probability of heads is $p$ and probability of tails is $q$).

The terms in the power represent the probabilities of getting heads or tails. And by writing out the product as a polynomial of $x$ we get the cases individual cases of the number of heads and tails.

Let's use as example n=3.

$$\begin{array}{rcl} (q+px)^3& =& (q+px)(q+px)(q+px) \\ &=& \overbrace{qqq + qqpx + qpxq + pxqq + qpxpx + pxqpx + pxpxq + pxpxpx}^{\substack{\text{8 terms for the 8 possible outcomes of heads and tails}\\\text{the number of $x$'s relate to the number of heads}}} \\ &=& \overbrace{(1 q^3) x^0 + (3 q^2p) x^1 + (3 q^1p^2) x^2 + (1 p^3) x^3 }^{\substack{\text{expression regrouped into groups with equal powers of $x$'s}}}\\ &=& \sum_{k=0}^3 a_k x^k = \sum_{k=0}^3 P(K=k) x^k \end{array}$$

It are these coefficients $a_k$ in the polynomial expansion that are of interest as they relate to the probability of $k$ heads. Computing the polynomial is the same as writing down all the possible combinations of heads and tails. You can see the power terms $x^k$ as keeping track of how many times you had heads.

---

Now if you would do a subtraction of two binomial terms, for instance with $m=n=3$ then you could use a product like

$$\begin{array}{rcl}(q_1+p_1x)^3(q_2+p_2x^{-1})^3& =& (q_1+p_1x)(q_1+p_1x)(q_1+p_1x)(q_2+p_2x^{-1})(q_2+p_2x^{-1})(q_2+p_2x^{-1}) \\ &=&q_1q_1q_1q_2q_2q_2 + q_1q_1q_1q_2q_2p_2/x + \text{ 62 other terms} \\&=& q_1^3p_2^3 x^{-3} + (3p_1q_1^2q_3^2+3 q_1^3p_2^2q_2^1) x^{-2} + (9p_1^1q_1^2p_2^2q_2^1+3 p_1^2q_1^1p_2^3 + 3 q_1^3p_2^1q_2^2) x^{-1} + (9 p_1^2q_1^1p_2^2q_2^1 + 9 p_1^1q_1^2p_2^1q_2^2 + p_1^3p_1^3 + q_1^3q_2^3) x^0 + (9p_1^2q_1^1p_2^1q_2^2+3 p_1^1q_1^2q_2^3 + 3 p_1^3p_2^2q_2^1) x^{1} + (3p_1^2q_1q_2^3+3 p_1^3p_2^1q_2^2) x^{2} + p_1^3q_2^3 x^{3} \\ &=& \sum_{k=-3}^3 a_k x^k = \sum_{k=-3}^3 P(K_1-K_2=k) x^k \end{array}$$

Now the power in the polynomial (or Laurent series since there are negative powers) represents the writing out of the product of the 6 terms as a sum of combinations of the 64 terms like $q_1q_1q_1q_2q_2p_2/x$ and keeps track how often in those terms we had a factor $x^{1}$ (relating to adding a heads) or a factor $x^{-1}$ (relating to subtracting a heads).

---

The coefficients for the series $a_k$ relate to the probabilities and are the values that you want to know. The formula with the derivatives are just a way to 'read out' those coefficients. The formula works when you do not have negative powers in the series. But you can solve this, as Whuber already commented, by multiplying the series with $x^m$.

---

The use of the probability generating function for this difference of two binomial distributions works, but it is not exciting. It becomes more interesting when we can perform some manipulations that simplify the function or expressions. An example is in this dice problem that can be solved by using a generating function: https://stats.stackexchange.com/a/492027/ where two tricks are applied (one is to rewrite a sum as a simple fraction $P(x) = \sum_{k=0}^\infty (5x+5x^2+5x^3+5x^4)^k = \frac{1-x}{1-6x+5x^5}$, a second is finding a recursive relation for the coefficients of the power series representation)