🌍从掷地球试验中学习贝叶斯推断

我们知道地球大约70%的面积被水覆盖，假设现在地球🌍像一颗小球一样在你手中，你将它向天空一抛，触地的那个点是水域的概率是多少？

本文资料和代码来自这里。

我们先假设我们完全不知道地球面积有多少的比例被水覆盖，但抛掷15次我们发现有8次是水，我们如何知道该结果的后验概率分布呢？

这里我们使用网格逼近，将0到1区间分割为1000个小区间，利用贝叶斯公式计算后验概率分布：

p_grid = seq(from=0, to=1, length.out=1000)
prior = rep(1, 1000)  # 无先验知识时，直接用均匀分布，这与p_grid一一对应
prob_data = dbinom(8, size=15, prob=p_grid)
posterior = prob_data * prior
posterior = posterior / sum(posterior)

然后我们从结果中抽样：

set.seed(100)
samples = sample(p_grid, prob=posterior, size=1e4, replace=TRUE)
mean(samples)

抽样的均值是0.53左右：

> mean(samples)
[1] 0.5281364

当我们已知地球大部分是水这一知识时结果会不会改善呢？我们修改先验知识重新进行计算：

p_grid = seq(from=0, to=1, length.out=1000)
prior = c(rep(0, 500), rep(1, 500))  # 已知水占比不可能低于0.5，所以直接给0值
prob_data = dbinom(8, size=15, prob=p_grid)
posterior = prob_data * prior
posterior = posterior / sum(posterior)

set.seed(100)
samples2 = sample(p_grid, prob=posterior, size=1e4, replace=TRUE)
mean(samples2)

结果：

> mean(samples2)
[1] 0.6054694

事实上，水占比为70%，我们将两次试验的结果以及真实值都在同一个图中画出来。

library(rethinking)

dens( samples , xlab="p" , xlim=c(0,1) , ylim=c(0,6))
dens( samples2 , add=TRUE, lty=2)
abline( v=0.7 , col="red")

现在有另外一个问题：我们想要精确地估计地球水域面积的占比，使得p值后验概率分布99%区间的区间都在0.05范围内，也就是上边界与下边界（99%区间）的差值是0.05，我们需要抛掷多少次地球？

下面创建一个函数用于计算给定抛掷次数p值后验分布的区间差值。

f = function(N) {
  p_true = 0.7
  W = rbinom(1, size = N, prob = p_true)
  p_grid = seq(from = 0, to = 1, length.out = 1000)
  prior = rep(1, 1000)
  prob_data = dbinom(W, size= N, prob = p_grid)
  posterior = prob_data * prior
  posterior = posterior / sum(posterior)
  samples = sample(p_grid, prob = posterior, size = 1e4, replace = TRUE)
  PI99 = PI(samples, 0.99)
  as.numeric(PI99[2] - PI99[1])
}

为了更好地对比，我们设置不同的抛掷次数，同一抛掷次数我们也计算多次，以观察计算结果的稳定性。

Nlist = c(20, 50, 100, 200, 500, 1000, 2000)
Nlist = rep(Nlist, each = 100)
width = sapply(Nlist, f)
plot(Nlist, width)
abline(h=0.05, col="red")

可以看到，大约需要2000次，我们可以将99%的区间差稳定在0.05左右。估计手很累吧？😋