机器算法验证 - 分层 Gamma-Poisson CDF？ - 吾爱随笔录

分层 Gamma-Poisson CDF？

机器算法验证贝叶斯泊松分布计算统计伽马分布

2022-04-03 09:15:41

评估 CDF 的计算效率最高的方法是什么

P (X \leq x ∣ r, v)

$P(X \leq x \mid r,v)$ 其中和

X \sim Poisson (λ)

$X \sim \operatorname{Poisson}(\lambda)$

λ \sim Gamma (r, v)

$\lambda \sim \operatorname{Gamma}(r,v)$

之后我看不到下一个明显的步骤

P (X \leq x ∣ r, v) = \int_{λ} P (X \leq x ∣ λ) P (λ ∣ r, v) d λ

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

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

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

我会很高兴有一个分析表达式，一个指向计算评估这个的最快方法的指针，或者最好的所有 Python 代码。

2个回答

如果你用另一种方式积分，你会得到一个封闭形式的表达式：

\begin{aligned} P (X \leq x | r, ν) & = \int_{0}^{\infty} P (X \leq x | λ) P (λ | r, v) d λ \\ = \int_{0}^{\infty} \sum_{k = 0}^{x} \frac{λ^{k} \exp (- λ)}{k!} \frac{λ^{r - 1} e^{- λ / v}}{Γ (r) v^{r}} d λ \\ = \frac{1}{Γ (r) v^{r}} \sum_{k = 0}^{x} \frac{1}{k!} \int_{0}^{\infty} λ^{k + r - 1} e^{- λ {1 + v^{- 1}}} d λ \\ = \sum_{k = 0}^{x} \frac{1}{k!} \frac{v^{- r}}{(1 + v^{- 1})^{k + r}} \frac{Γ (k + r)}{Γ (r)} \end{aligned}

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

是的，当没有 PDF 的封闭形式时，您可以使用 Metropolis-Hastings 算法。一旦你理解了这个算法，你就可以很容易地自己编程，但是在谷歌搜索中，很多人已经在 python 中这样做了。这里 http://www.nehalemlabs.net/prototype/blog/2014/02/24/an-introduction-to-the-metropolis-method-with-python/。

我在 R 中使用以下代码和一些添加的图完成了此操作：

## 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)

R 代码很笨拙，但出于说明目的而编写。

其它你可能感兴趣的问题

上一篇使用混合效应模型的真正动机，以及何时使用它们以及何时不使用它们下一篇测试趋势（序数预测变量，连续结果）？