python - 如何建立三次抛硬币的贝叶斯仿真模型

Question

假设我们抛了 8 次有偏差的硬币(我们不知道它有多大偏差)，到目前为止我们记录了 5 次正面(H)到 3 次反面(T)。接下来的 3 次抛掷都是反面的概率是多少？换句话说，我们想知道在第 11 次 throw 后出现 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 分布不同。