# How can I apply the Poisson ($\mu$) distribution to two series of random draws?

I have the following question:

A box contains 1000 balls, of which 2 are black and the rest are white. If two series of 1000 draws are made at random from this box, what approximately, is the chance that they produce the same number of black balls?

I believe it an exercise in applying The Poisson ($$\mu$$) Distribution: The Poisson distribution with parameter $$\mu$$ is the distribution of probabilities $$P_{\mu}(k)$$ over $$\{0,1,2,...\}$$ defined by $$P_{\mu}(k)=e^{-\mu}\frac{\mu^k}{k!}$$

Let $$\mu = 1000 \cdot \frac{2}{1000} = 2$$ so that $$P_{\mu}(k)=\sum_{k=0}^{\infty}P(k)\cdot P(k)$$ $$P_{\mu}(k)= \sum_{k=0}^{\infty} e^{-2}\frac{2^k}{k!} \cdot e^{-2}\frac{2^k}{k!}$$ $$=e^{-4}\cdot \sum_{k=0}^{\infty}(\frac{2^k}{k!})^2$$

And then I look up the answer and see that $$\sum_{k=0}^{\infty}(\frac{2^k}{k!})^2$$ translates to $$I_0(4)$$ which is a modified Bessel function. However I havn't learnt this, and Bessel isn't in the index of my textbook.

I thought maybe $$\sum_{k=0}^{\infty}(\frac{2^k}{k!})^2 = (e^2)^2$$ but that would make the whole thing =1... which maybe makes sense since they question didn't specify with replacement.

• You can't interchange summation and squaring. – Glen_b Oct 16 at 4:26

Comments: I'm not sure why you think you're supposed to use a Poisson approximation.

Poisson approximation. With your Poisson method, approximating the (infinite?) sum with R:

f = dpois(0:10000, 2)
sum(f^2)
[1] 0.2070019


But I'm not sure why you'd sum the series beyond 1000 draws?

Actually, almost all of the probability is in the first few terms:

f = dpois(0:10, 2);  f.2 = f^2;  cf.2 = cumsum(f.2)
cbind(f, f.2, cf.2)  # 'binds' column vectors

f          f.2       cf.2
[1,] 1.353353e-01 1.831564e-02 0.01831564
[2,] 2.706706e-01 7.326256e-02 0.09157819
[3,] 2.706706e-01 7.326256e-02 0.16484075
[4,] 1.804470e-01 3.256114e-02 0.19740189
[5,] 9.022352e-02 8.140284e-03 0.20554217
[6,] 3.608941e-02 1.302445e-03 0.20684462
[7,] 1.202980e-02 1.447162e-04 0.20698933
[8,] 3.437087e-03 1.181356e-05 0.20700114
[9,] 8.592716e-04 7.383478e-07 0.20700188
[10,] 1.909493e-04 3.646162e-08 0.20700192
[11,] 3.818985e-05 1.458465e-09 0.20700192


Simulating with 10 million iterations, result is $$0.2071 \pm 0.0013.$$

set.seed(1016)
x = rpois(10^7, 2)
y = rpois(10^7, 2)
mean(x==y)
[1] 0.207069      # aprx P(Same twice) = 0.2070019
2*sd(x--y)/sqrt(10^7)
[1] 0.001264853   # aprx 95% margin of simulation error


So provided the Poisson model is correct, the answer to four places is 0.2070.

Hypergeometric. If you choose without replacement, then you'll always get exactly 2 black balls.

Binomial. If you choose with replacement, the exact model is binomial.

f = dbinom(0:1000, 1000, 2/1000)
sum(f^2)
[1] 0.2071337


Maybe there is a nice binomial identity for summing the 1001 terms.

Simulation of binomial model gives $$P(\text{Both Same}) = 0.2073 \pm 0.0003,$$ which includes the exact value.

set.seed(1017)
x = rbinom(10^7, 1000, 2/1000)
y = rbinom(10^7, 1000, 2/1000)
mean(x==y)
[1] 0.207295
2*sd(x==y)/sqrt(10^7)
[1] 0.0002563777

• This is the correct answer. For the first method (Poisson Approximation) - how could this calculation be done with a scientific calculator? Or would there be a table that would need to be consulted? Previous to this question I've figured out the probabilities of 1) Less than 2 black ball, 2) exactly 2 black balls, 3) more than 2 black balls - if that helps. – Elliott de Launay Oct 16 at 13:57
• On a calculator, if you had the patience to sum only 1st 11 terms of series, you'd get close: In R, f = dpois(0:10, 2); sum(f^2) returns 0.2070019 because the first few terms are the only large ones. You could see when the squares get too small to matter. // But you should start learning R if this isn't your last prob/stat course. Not all of R. Maybe just a few dozen functions at the start, mostly intuitive. – BruceET Oct 16 at 19:23
• FWIW, the Binomial answer can be expressed (with $n=1000$ and $p=2/1000$) as $(1-p)^{2n}\,_2F_1(-n;-n;1;p^2/(1-p)^2).$ It expands to a polynomial of degree $2n.$ Indeed, using the hypergeo library R (a "scientific calculator" if there ever was one!) will compute n <- 1000; p <- 2/1000; ((1-p)^(2*n) * hypergeo(-n,-n,1,(p/(1-p))^2)) as 0.2071337, agreeing with the dbinom value to 15 significant digits. Unfortunately there's only about 10% improvement in computation time. – whuber Oct 16 at 20:13
• FWW (again, now that we're edging toward giving final hwk solutions): For the first few terms of binomial: f = dbinom(0:10, 1000, 2/1000); sum(f^2) returns 0.2071337. // IMHO unless undisclosed directions mandate Poisson aprx, I see no reason to avoid exact binomial computation. – BruceET Oct 16 at 21:52