Compared Two Poisson Variables

Background

Last Friday, I had to submit a homework which required me to evaluate $\Pr(A > B)$ and $\Pr(A = B)$, where $A$ and $B$ were two independent Poisson random variables with parameters $\alpha$ and $\beta$ respectively.

Problem

I then started evaluating the sum.

\[ \Pr(A > B) = \sum_{i = 1}^\infty \sum_{j = 0}^{i - 1} \frac{e^{-\alpha} \alpha^i}{i!} \cdot \frac{e^{-\beta} \beta^j}{j!} \]

Then I was stuck. I couldn’t compute this sum also.

\[ \Pr(A = B) = \sum_{i = 0}^\infty \frac{e^{-(\alpha + \beta)} \alpha^i \beta^i}{(i!)^2} \]

Fact

I googled for a solution for hours, and after I saw equation (3.1) in a paper, I gave up finding exact solutions.¹ As a supporter of free software, I avoided using M$ Ex*, and wrote a program in C++ to approximate the above probabitities by directly adding them term by term.

Source code

#include <iostream> #include <cmath> using namespace std; double pXy(double x, double y, int N); double pxy(double x, double y, int N); int main(void) { double pAb,paB,pab,a,b; int N; cout << "Assume that Poisson r.v. A and B are indepedent" << endl; cout << "Parameter for A: "; cin >> a; cout << "Parameter for B: "; cin >> b; cout << "Number of terms to be added (100 <= N <= 1000): "; cin >> N; pAb = pXy(a,b,N); paB = pXy(b,a,N); pab = pxy(a,b,N); cout << "P(A > B) = " << pAb << ", P(A < B) = " << paB << ", P(A = B) = " << pab << endl; } /* P(X > Y) */ double pXy(double x, double y, int N) { double ans = 0; for (int i = 1; i <= N; i++) { for (int j = 0; j < i ; j++) { double term = 1; for (int k = 1; k <= i; k++) term *= x / k; for (int k = 1; k <= j; k++) term *= y / k; ans += term; } } return ans * exp(-x - y); } /* P(X = Y) */ double pxy(double x, double y, int N) { double ans = 0; for (int i = 0; i <= N; i++) { double term = 1; for (int k = 1; k <= i; k++) term *= x / k * y / k; ans += term; } return ans * exp(-x - y); }

Sample output

Assume that Poisson r.v. A and B are indepedent
Parameter for A: 1.6
Parameter for B: 1.4
Number of terms to be added (100 <= N <= 1000): 8
P(A > B) = 0.423023, P(A < B) = 0.335224, P(A = B) = 0.241691

Lessons learnt

A one-line method for writing the content of a function which returns the factorial of a number.

URL: http://progopedia.com/example/factorial/
Evaluation of a function inside GDB

URL: http://stackoverflow.com/q/1354731/

Keller, J. B. (1994). A characterization of the Poisson distribution and the probability of winning a game. The American Statistician, 48(4), 294–298. ↩

Blog 1

Random Talk on Random Thoughts