这是关于 Jim Albert 的“使用 R 的贝叶斯计算”中的一个练习的初学者问题。请注意，虽然这可能是家庭作业，但在我的情况下并非如此，因为我正在学习 R 中的贝叶斯方法，因为我认为我可能会在未来的分析中使用它。

无论如何，虽然这是一个具体的问题，但它可能涉及对贝叶斯方法的基本理解。

因此，在练习 2.2 中，Jim Albert 要求我们分析投掷硬币的实验。看这里。我们将使用先验直方图，即将可能p值的空间划分为 10 个长度间隔，.1并为这些间隔分配先验概率。

由于我知道真正的概率是.5，并且我认为宇宙改变概率定律或一分钱坚固耐用的可能性很小，我的先验是：

prior <- c(1,5,20,100,5000,5000,100,20,5,1)
prior <- prior/sum(prior)

沿区间中点

midpt <- seq(0.05, 0.95, by=0.1)

到目前为止，一切都很好。接下来，我们将硬币旋转 20 次并记录成功（正面）和失败（尾部）的次数。轻松完成：

y <- rbinom(n=20,p=.5,size=1)
s <- sum(y==1)
f <- sum(y==0)

在我的实验中，s == 7并且f == 13. 接下来是我不明白的部分：

通过 (1) 在 (0,1) 上的值的网格上计算 p 的后验密度和 (2) 从网格中获取具有替换的模拟样本，从后验分布进行模拟。（函数 histprior和sample有助于此计算）。根据您的数据，区间概率如何变化？

这是这样做的：

p <- seq(0,1, length=500)
post <- histprior(p,midpt,prior) * dbeta(p,s+1,f+1)
post <- post/sum(post)

ps <- sample(p, replace=TRUE, prob = post)

但是我们为什么要这样做呢？

我们可以很容易地通过将先验乘以适当的可能性来获得后验密度，就像上面块的第二行所做的那样。这是后验分布图：

由于后验分布是有序的，我们可以通过汇总后验密度的元素来获得直方图先验中引入的区间的结果：

post.vector <- vector()
post.vector[1] <- sum(post[p < 0.1])
post.vector[2] <- sum(post[p > 0.1 & p <= 0.2])
post.vector[3] <- sum(post[p > 0.2 & p <= 0.3])
post.vector[4] <- sum(post[p > 0.3 & p <= 0.4])
post.vector[5] <- sum(post[p > 0.4 & p <= 0.5])
post.vector[6] <- sum(post[p > 0.5 & p <= 0.6])
post.vector[7] <- sum(post[p > 0.6 & p <= 0.7])
post.vector[8] <- sum(post[p > 0.7 & p <= 0.8])
post.vector[9] <- sum(post[p > 0.8 & p <= 0.9])
post.vector[10] <- sum(post[p > 0.9 & p <= 1])

（R 专家可能会找到一种更好的方法来创建该向量。我想这可能与sweep？）

round(cbind(midpt,prior,post.vector),3)

      midpt prior post.vector
 [1,]  0.05 0.000       0.000
 [2,]  0.15 0.000       0.000
 [3,]  0.25 0.002       0.003
 [4,]  0.35 0.010       0.022
 [5,]  0.45 0.488       0.737
 [6,]  0.55 0.488       0.238
 [7,]  0.65 0.010       0.001
 [8,]  0.75 0.002       0.000
 [9,]  0.85 0.000       0.000
[10,]  0.95 0.000       0.000

此外，我们有 500 个来自后验分布的图，它们告诉我们没有什么不同。这是模拟平局的密度图：

现在我们使用模拟数据通过计算区间内模拟的比例来获得区间的概率：

sim.vector <- vector()
sim.vector[1] <- length(ps[ps < 0.1])/length(ps)
sim.vector[2] <- length(ps[ps > 0.1 & ps <= 0.2])/length(ps)
sim.vector[3] <- length(ps[ps > 0.2 & ps <= 0.3])/length(ps)
sim.vector[4] <- length(ps[ps > 0.3 & ps <= 0.4])/length(ps)
sim.vector[5] <- length(ps[ps > 0.4 & ps <= 0.5])/length(ps)
sim.vector[6] <- length(ps[ps > 0.5 & ps <= 0.6])/length(ps)
sim.vector[7] <- length(ps[ps > 0.6 & ps <= 0.7])/length(ps)
sim.vector[8] <- length(ps[ps > 0.7 & ps <= 0.8])/length(ps)
sim.vector[9] <- length(ps[ps > 0.8 & ps <= 0.9])/length(ps)
sim.vector[10] <- length(ps[ps > 0.9 & ps <= 1])/length(ps)

（再次：有没有更有效的方法来做到这一点？）

总结结果：

round(cbind(midpt,prior,post.vector,sim.vector),3)

      midpt prior post.vector sim.vector
 [1,]  0.05 0.000       0.000      0.000
 [2,]  0.15 0.000       0.000      0.000
 [3,]  0.25 0.002       0.003      0.000
 [4,]  0.35 0.010       0.022      0.026
 [5,]  0.45 0.488       0.737      0.738
 [6,]  0.55 0.488       0.238      0.236
 [7,]  0.65 0.010       0.001      0.000
 [8,]  0.75 0.002       0.000      0.000
 [9,]  0.85 0.000       0.000      0.000
[10,]  0.95 0.000       0.000      0.000

毫无疑问，simmultion 除了基于它的后验之外没有产生其他结果。因此，为什么我们首先要绘制这些模拟？