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

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).

• The PGF is a formal Laurent series. You don't need to think of it as a function of $z.$ If you really must, then simply multiply $G_S$ by $z^n$ and analyze that.
– whuber
Commented Apr 21, 2023 at 15:12
• Ok thanks! I'll try that. Commented Apr 21, 2023 at 15:31

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

• Could you explain how this answers the question?
– whuber
Commented Apr 21, 2023 at 15:12
• Whuber, check the edited post. I guess this is probably in the same vein as to what you remarked above. I found the reference helpful but that was marred with typos. Commented Apr 21, 2023 at 19:40
• It's unclear which sequence you're discussing. As I remarked in another comment, you don't need any Real Analysis for this, because these series can be treated purely formally: the rules of multiplying these formal power series are equivalent to the rules for summing random variables supported on the integers and the derivative can be taken purely formally, too.
– whuber
Commented Apr 21, 2023 at 20:18
• I agree completely. As for the sequence, they were the values taken by $X$ in an increasing order. Commented Apr 21, 2023 at 20:39
• I see, thank you. But there's no need to take limits. Without any loss of generality you can analyze the shifted variable $X - \alpha_{-n},$ whose smallest value is almost surely $0.$ In the case of the question, where all values are integral, you thereby arrive at a polynomial and the question is a purely algebraic one.
– whuber
Commented Apr 21, 2023 at 20:49

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)