R-基本分层贝叶斯分析的代码

Question

我有此代码写在winbugs中：

n <- 100
x1 <- rbinom(n,1,.7)
x2 <- rbinom(n,1,.5)
sum(x1)
sum(x2)
model{
x1 ~ dbin(p1, n) x2 ~ dbin(p2, n) p1 ~ dbeta(a1, b1) p2 ~ dbeta(a2,b2)
diff <- p1 - p2 p.value <- step(diff)
} list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)

我在如何在R/JAGS中执行此操作时遇到了麻烦。实际上，我什至不完全确定该代码要做什么（我认为计算后代？）。我以前从未使用过winbugs，也是R的新手。这也是我的第一个贝叶斯课，一旦引入了代码，我就会迷路。

此外，我如何计算比例差异的后均值和标准偏差？或如何找到p1大于p2的后验概率，以及是否重要？

Answer 1

如果您需要帮助 rjags 设置，我希望您的TA或同学可以为您提供帮助。该答案重点是解释代码的作用。

出于教学目的，让我们对此数据强加一个简单的叙述。假设我们有两个硬币，我们想知道一个硬币是否比另一个更有可能翻转头。代码分为两个部分。

生成/正向模拟

第一个是一个正向模拟，我们已经知道两枚硬币翻转头的概率是什么。

## how many times we will flip each coin
n <- 100
## coin 1 has 70% chance to land on heads
## simulate n flips
x1 <- rbinom(n, 1, .7)
## coin 2 has 50% chance to land on heads
## simulate n flips
x2 <- rbinom(n, 1, .5)
## count number of coin 1 heads
sum(x1) # 70
## count number of coin 2 heads
sum(x2) # 54

后验概率

现在，我们采用这些模拟数据并尝试逆转实验。也就是说，如果我们从一个硬币中观察到70个头，另一个来自另一个硬币的头，我们可以说出每个硬币必须翻转头部的概率吗？具体而言，代码将要问："在一个统一的假设下，关于每个硬币具有的真实概率的均匀假设，硬币1比硬币2更可能翻转头部的（后验）概率是什么？"

计算似然

jags这样做的方式是对MCMC，其中样品是从所有可能的参数配置的空间（在这种情况下，是p1和p2的值），加权每个配置生成观察到的数据的可能性。因此，jags代码需要完成的主要事情是定义如何生成参数值，以及如何计算这些值，如何计算数据的可能性。这是模型代码的第一部分。

兴趣的计算值

模型部分中代码的第二部分正在测试我们提出的问题。这是通过计算一些不影响配置可能性的其他变量来完成的。相反，他们告诉我们有关要采样的配置的信息。具体而言，diff将跟踪两个硬币概率之间的差异分布，并且p.value将跟踪p1大于p2的频率。

model{
## COMPUTE LIKELIHOODS
## compute likelihood that heads resulted from coins
## with given probabilities after `n` coin flips
x1 ~ dbin(p1, n) 
x2 ~ dbin(p2, n) 
## SAMPLE PARAMETERS
## randomly select coin probabilities from Beta distributions
## Note that since these are all 1, it's really just a Uniform[0,1]
p1 ~ dbeta(a1, b1)
p2 ~ dbeta(a2, b2)
## HYPOTHESIS TEST
## compute how often coin1's probability is greater than coin2's
diff <- p1 - p2 
p.value <- step(diff)
} 
## INPUT VALUES
list(n = 100, x1 = 70, x2 = 54, a1 = 1, b1 = 1, a2 = 1, b2 = 1)