5
$\begingroup$

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.

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

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*}

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

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.