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.
 A: 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*}
A: 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)
    ## deciding to jump to lambda cand or not
    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. 
