## parameters
set.seed(1)
ns <- 10^5
nbinom <- 40
pbinom <- 0.5

### generate sample
Y <- rbinom(ns, nbinom, pbinom)
X <- rpois(ns, Y)

### plot histogram
h <- hist(X, breaks = seq(0,max(X)+1)-0.5, freq = 0)

### add curve
xs <- seq(0,max(X)) 
#lines(xs,dpois(xs,nbinom*pbinom)*ns, col = 2)
lines(xs,dnorm(xs+0.5, mean = nbinom*pbinom, 
                     sd = sqrt(nbinom*pbinom*(2-pbinom))), col = 2)

### binomial approximation
mu = nbinom*pbinom
vr = (nbinom*pbinom*(2-pbinom))
p = 1-mu/vr
r = mu*(1-p)/p
lines(xs,dnbinom(xs, size = r, p= 1-p ), col = 3)

legend(max(X),max(h$density), c("normal approximation", "negative binomial approximation"),
       cex= 0.7, col = c(2,3),lty = 1, xjust = 1)

## parameters
set.seed(1)
ns <- 10^5
nbinom <- 40
pbinom <- 0.5

### generate sample
Y <- rbinom(ns, nbinom, pbinom)
X <- rpois(ns, Y)

### plot histogram
h <- hist(X, breaks = seq(0,max(X)+1)-0.5, freq = 0)

### add curve
xs <- seq(0,max(X)) 
#lines(xs,dpois(xs,nbinom*pbinom)*ns, col = 2)
lines(xs,dnorm(xs+0.5, mean = nbinom*pbinom, 
                     sd = sqrt(nbinom*pbinom*(2-pbinom))), col = 2)

### binomial approximation
mu = nbinom*pbinom
vr = (nbinom*pbinom*(2-pbinom))
p = 1-mu/vr
r = mu*(1-p)/p
lines(xs,dnbinom(xs, size = r, p= 1-p ), col = 3)

legend(max(X),max(h$density), c("normal approximation", "negative binomial approximation"),
       cex= 0.7, col = c(2,3),lty = 1, xjust = 1)

Normal approximation

One approach is to use a normal approximation by matching the moments.

The compound Poisson distribution is a mixture distribution of Poisson distributions with weights distributed according to a binomial distribution $w_i = P_{binom}(i\vert n,p)$

$$P(x) = \sum_{i=0}^n w_iP_{Poisson}(x\vert \lambda = i )$$

• The mean of the distribution will be equal to $$E[X] = \sum_{i=0}^n w_i E[X\vert i] = \sum w_i i = np$$ the mean of the binomial distribution.
• The variance will be $$\begin{array}{} E[(X-\mu_X)^2] &=& \sum_{i=0}^n w_i (\sigma_i^2 + \mu_i^2 - \mu^2) \\& = &\sum_{i=0}^n w_i (i + i^2 - (np)^2) \\&=& \sum_{i=0}^n w_i (i + i^2) - \sum_{i=0}^n (np)^2\\ &=& np(1-p) + (np) \end{array}$$ the variance plus the mean of the binomial distribution

Computational

The example below is computed with $n=40$ and $p=0.5$ giving $\mu_X = 20$ and $\sigma_X^2 = 30$

### parameters
set.seed(1)
ns <- 10^5
nbinom <- 40
pbinom <- 0.5

### generate sample
X <- rbinom(ns, nbinom, pbinom)
Y <- rpois(ns, X)

### plot histogram
hist(Y, breaks = seq(0,max(Y)+1)-0.5, freq = 0)

### add curve
xs <- seq(0,max(Y)) 
#lines(xs,dpois(xs,nbinom*pbinom)*ns, col = 2)
lines(xs,dnorm(xs+0.5, mean = nbinom*pbinom, 
                     sd = sqrt(nbinom*pbinom*(2-pbinom))), col = 2)

The problem with the normal approximation is that it does not work so well when $\sigma_X$ is not a lot smaller than $\mu_X$. This occurs when $\mu_X = np$ is small.

See for instance the plot below with $n=4000$ and $p=0.0005$ giving $\mu_X = 2$ and $\sigma_X^2 = 1.99975$

Normal approximation

One approach is to use a normal approximation by matching the moments.

The compound Poisson distribution is a mixture distribution of Poisson distributions with weights distributed according to a binomial distribution $w_i = P_{binom}(i\vert n,p)$

$$P(x) = \sum_{i=0}^n w_iP_{Poisson}(x\vert \lambda = i )$$

• The mean of the distribution will be equal to $$E[X] = \sum_{i=0}^n w_i E[X\vert i] = \sum w_i i = np$$ the mean of the binomial distribution.
• The variance will be $$\begin{array}{} E[(X-\mu_X)^2] &=& \sum_{i=0}^n w_i (\sigma_i^2 + \mu_i^2 - \mu^2) \\& = &\sum_{i=0}^n w_i (i + i^2 - (np)^2) \\&=& \sum_{i=0}^n w_i (i + i^2) - \sum_{i=0}^n (np)^2\\ &=& np(1-p) + (np) \end{array}$$ the variance plus the mean of the binomial distribution

Negative binomial approximation

We can do the same sort of approximation (matching the moments) with a negative binomial distribution.

In the examples below, we have plotted this approximation as well. The curve resembles more closely the histogram than the normal approximation, but it is still not exactly the same.

Computational

The example below is computed with $n=40$ and $p=0.5$ giving $\mu_X = 20$ and $\sigma_X^2 = 30$

## parameters
set.seed(1)
ns <- 10^5
nbinom <- 40
pbinom <- 0.5

### generate sample
Y <- rbinom(ns, nbinom, pbinom)
X <- rpois(ns, Y)

### plot histogram
h <- hist(X, breaks = seq(0,max(X)+1)-0.5, freq = 0)

### add curve
xs <- seq(0,max(X)) 
#lines(xs,dpois(xs,nbinom*pbinom)*ns, col = 2)
lines(xs,dnorm(xs+0.5, mean = nbinom*pbinom, 
                     sd = sqrt(nbinom*pbinom*(2-pbinom))), col = 2)

### binomial approximation
mu = nbinom*pbinom
vr = (nbinom*pbinom*(2-pbinom))
p = 1-mu/vr
r = mu*(1-p)/p
lines(xs,dnbinom(xs, size = r, p= 1-p ), col = 3)

legend(max(X),max(h$density), c("normal approximation", "negative binomial approximation"),
       cex= 0.7, col = c(2,3),lty = 1, xjust = 1)

The problem with the normal approximation is that it does not work so well when $\sigma_X$ is not a lot smaller than $\mu_X$. This occurs when $\mu_X = np$ is small.

See for instance the plot below with $n=4000$ and $p=0.0005$ giving $\mu_X = 2$ and $\sigma_X^2 = 1.99975$