# Hierarchical Gamma-Poisson CDF?

What is the most computationally efficient way to evaluate the CDF

$$P(X \leq x \mid r,v)$$ where $$X \sim \operatorname{Poisson}(\lambda)$$ and $$\lambda \sim \operatorname{Gamma}(r,v)$$

I can't see the next obvious step after

$$P(X \leq x\mid r,v) = \int_\lambda P(X \leq x\mid \lambda)P(\lambda\mid r,v) \,d\lambda$$

$$=\int_\lambda \frac{\Gamma(x+1,\lambda)}{x!} \frac{\lambda^{r-1}e^{-\lambda/v}}{\Gamma(r)v^r} \, d\lambda$$

I would be happy with either an analytical expression, a pointer for the fastest way to computationally evaluate this or best of all Python code that already does.

• It should give you a negative binomial, shouldn't it? – Glen_b -Reinstate Monica Jan 8 '15 at 3:42

If you integrate the other way, you get a closed form expression: \begin{align*} P(X \leq x|r,\nu) &= \int_0^\infty P(X \leq x| \lambda)P(\lambda|r,v) \text{d}\lambda\\ &= \int_0^\infty \sum_{k=0}^x \frac{\lambda^k\exp(-\lambda)}{k!} \frac{\lambda^{r-1}e^{-\lambda/v}}{\Gamma(r)v^r} \text{d}\lambda\\ &= \frac{1}{\Gamma(r)v^r} \sum_{k=0}^x \frac{1}{k!} \int_0^\infty \lambda^{k+r-1}e^{-\lambda\{1+v^{-1}\}}\text{d}\lambda\\ &= \sum_{k=0}^x \frac{1}{k!} \frac{v^{-r}}{(1+v^{-1})^{k+r}}\,\frac{\Gamma(k+r)}{\Gamma(r)} \end{align*}

• Thanks, I had a feeling there was a closed form but I couldn't see how to get there! – rwolst Jan 8 '15 at 17:40

Yes when there is no closed form of the PDF you can use the Metropolis-Hastings Algorithm. Once you understand the algorithm you can easily program it yourself but upon a google search many have done so in python. Here http://www.nehalemlabs.net/prototype/blog/2014/02/24/an-introduction-to-the-metropolis-method-with-python/.

I have done this in R with the following code and some added plots:

## M-H algorithm to draw values from a Poisson with lambda=8
## using a negative binomial distribution with mean equal to
## lambda current as the proposal dist.
## intial values and vectors to start the loop, algorithm uses
## three chains
nsim <- 1000
lam.vec1 <- numeric(nsim)
lam.vec2 <- numeric(nsim)
lam.vec3 <- numeric(nsim)
lam.vec1[1] <- 7
lam.vec2[1] <- 6
lam.vec3[1] <- 5
jump.vec1 <- numeric(nsim - 1)
jump.vec2 <- numeric(nsim - 1)
jump.vec3 <- numeric(nsim - 1)
s <- 4
for (i in 2:nsim) {
## assigning lambda current for the three chains
lam.curr1 <- lam.vec1[i - 1]
lam.curr2 <- lam.vec2[i - 1]
lam.curr3 <- lam.vec3[i - 1]
## obtaining lambda candidate for three chains by drawing from
## the proposal dist.
lam.cand1 <- rnbinom(1, size = s, mu = lam.curr1)
lam.cand2 <- rnbinom(1, size = s, mu = lam.curr2)
lam.cand3 <- rnbinom(1, size = s, mu = lam.curr3)
## numerator for the M-H ratio
r.num1 <- dpois(lam.cand1, lambda = 8) * dnbinom(lam.curr1,
size = s, mu = lam.cand1)
r.num2 <- dpois(lam.cand2, lambda = 8) * dnbinom(lam.curr2,
size = s, mu = lam.cand2)
r.num3 <- dpois(lam.cand3, lambda = 8) * dnbinom(lam.curr3,
size = s, mu = lam.cand3)
## denominator for M-H ratio
r.den1 <- dpois(lam.curr1, lambda = 8) * dnbinom(lam.cand1,
size = s, mu = lam.curr1)
r.den2 <- dpois(lam.curr2, lambda = 8) * dnbinom(lam.cand2,
size = s, mu = lam.curr2)
r.den3 <- dpois(lam.curr3, lambda = 8) * dnbinom(lam.cand3,
size = s, mu = lam.curr3)
## M-H ratio
r1 <- r.num1/r.den1
r2 <- r.num2/r.den2
r3 <- r.num3/r.den3
## accept with probability min(1,r1)
p.accept1 <- min(1, r1)
p.accept2 <- min(1, r2)
p.accept3 <- min(1, r3)
u.vec <- runif(3)
ifelse(u.vec[1] <= p.accept1, lam.vec1[i] <- lam.cand1, lam.vec1[i] <- lam.curr1)
ifelse(u.vec[2] <= p.accept2, lam.vec2[i] <- lam.cand2, lam.vec2[i] <- lam.curr2)
ifelse(u.vec[3] <= p.accept3, lam.vec3[i] <- lam.cand3, lam.vec3[i] <- lam.curr3)
## recording whether we jumped or not
jump.vec1[i - 1] <- ifelse(u.vec[1] <= p.accept1, 1, 0)
jump.vec2[i - 1] <- ifelse(u.vec[2] <= p.accept2, 1, 0)
jump.vec3[i - 1] <- ifelse(u.vec[3] <= p.accept3, 1, 0)
}
## Look at chains
plot(seq(1:nsim), lam.vec1, type = "l", ylab = expression(lambda),
col = 4, main = "Traceplot of Three Chains")
lines(seq(1:nsim), lam.vec2, col = 2)
lines(seq(1:nsim), lam.vec3, col = 3)
# mean(jump.vec1) mean(jump.vec2) mean(jump.vec3)
post.lam <- c(lam.vec1, lam.vec2, lam.vec3)
x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
plot(table(post.lam)/length(post.lam), col = 3, type = "h", xlab = expression(lambda),
ylab = "Density")
points(dpois(rep(1:20, 1), lambda = 8), lwd = 2)
legend("topright", legend = c("Analytic Poission"), pch = 1,
lwd = 2, col = c(1), bty = "n", cex = 1)


The R code is clunky but was written as such for illustrative purposes.