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

Question

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.

Answer 1

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