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

Question

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

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

最常见的解决方案：

如果我们考虑单次抛掷的概率： 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 分布不同。

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

How to build a Bayesian simulation model for flipping coin three times

python

bayesian

pymc3

mcmc

解析解