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

Question

想象一下，我们扔了一个有偏见的硬币8次（我们不知道它是多么有偏见），到目前为止我们记录了5个头（H）到3个尾（T）。接下来的3次投掷都是尾巴的概率是多少？换句话说，我们想知道在第11次投掷后有5Hs和6Ts的预期概率。

我想使用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分布不同。