如何建立抛硬币三次的贝叶斯模拟模型

Question

假设我们将一枚有偏差的硬币扔了 8 次（我们不知道它有多大偏差），到目前为止，我们记录了 5 个正面 (H) 到 3 个反面 (T)。接下来的 3 次投掷都是反面的概率是多少？换句话说，我们想知道在第 11 次投掷后出现 5H 和 6T 的预期概率。

我想使用 pyMC3 构建一个 MCMC 仿真模型来找到贝叶斯解决方案。针对这个问题，贝叶斯方法中也有一个解析解。因此，我将能够比较从模拟、分析方法以及经典最常见方法得出的结果。让我简要解释一下到目前为止我能做什么：

1. 最常见的解决方案：

如果我们考虑单次投掷的概率： E(T) = p = (3/8) = 0,375 那么，最终的答案是 E({T,T,T}) = p^3 = (3/8)^3 = 0,052。

2.1。贝叶斯解析解法：

请假设未知参数“p”代表硬币的偏差。 如果我们考虑单次投掷的概率： E(T) = Integral0-1[p * P(p | H = 5, T = 3) dp] = 0,400（我经过一些代数运算计算了结果） 同样，最终的答案是： E({T,T,T}) = Integral0-1[p^3 * P(p | H = 5, T = 3) dp] = 10/11 = 0,909。

2.2。使用 MCMC 模拟的贝叶斯解决方案： 当我们考虑单次投掷的概率时，我在pyMC3中建立了如下模型：

Head: 0 
Tail: 1
data = [0, 0, 0, 0, 0, 1, 1, 1]
import pymc3 as pm

with pm.Model() as coin_flipping:
    p = pm.Uniform(‘p’, lower=0, upper=1)
    y = pm.Bernoulli(‘y’, p=p, observed=data)
    trace = pm.sample(1000)
    pm.traceplot(trace)

运行此代码后，我得到后验平均值为 E(T) =0,398，这与解析解 (0,400) 的结果非常接近。到目前为止我很高兴，但这不是最终的答案。我需要一个模拟 E({T,T,T}) 概率的模型。如果有人在这一步上帮助我，我将不胜感激。

Answer 1

一种方法是根据经验使用 PyMC3 的后验预测采样。也就是说，一旦有了后验抽样，就可以从模型的随机参数化中生成抽样。 pymc3.sample_posterior_predictive() 方法将生成与原始观察数据大小相同的新样本。由于您只对三个翻转感兴趣，我们可以忽略它产生的额外翻转。例如，如果你想要 10000 组随机预测翻转，你会这样做

with pm.Model() as coin_flipping:
    # this is still uniform, but I always prefer Beta for proportions
    p = pm.Beta(‘p’, alpha=1, beta=1)

    pm.Bernoulli(‘y’, p=p, observed=data)

    # chains looked a bit waggly at 1K; 10K looks smoother
    trace = pm.sample(10000, random_seed=2019, chains=4)

    # note this generates (10000, 8) observations
    post_pred = pm.sample_posterior_predictive(trace, samples=10000, random_seed=2019)

然后看看接下来的三个翻转的频率是多少 (1,1,1)，我们可以这样做

np.mean(np.sum(post_pred['y'][:,:3], axis=1) == 3)
# 0.0919

解析解

在这个例子中，由于我们有一个分析后验预测分布（Beta-Binomial[k | n, a=4, b=6] - 详情请参阅the Wikipedia table of conjugate distributions），我们可以准确计算在接下来的三个翻转中观察到三个尾部的概率如下：

from scipy.special import comb, beta as beta_fn

n, k = 3, 3  # flips, tails
a, b = 4, 6  # 1 + observed tails, 1 + observed heads

comb(n, k) * beta_fn(n + a, n - k + b) / beta_fn(a, b)
# 0.09090909090909091

请注意，beta_fn 是 Euler Beta function，与 Beta 版不同。