发现贝叶斯的乐高积木

原文：https://towardsdatascience.com/https-medium-com-hankroark-finding-bayesian-legos-part1-b8aeb886afba

 

照片来源：FrédériqueVoisin-Demery / Flickr（CC BY 2.0）

我有一个很好的朋友Joe，本周路过他家时我顺便造访了他家。像平常一样，我们聊了天气（在太平洋西北地区已经比正常情况更热），新闻（主要是关于我们该如何应对炎热天气）和我们的孩子。我们俩都有孩子，而且他们都非常喜欢玩乐高积木。家里散落着大量的乐高积木就不可避免地会出现踩踏乐高积木的强烈痛苦，通常是在半夜或早上醒来制作咖啡时。为了避免出现那样的痛苦，我们就追在孩子的后面收拾被他们乱扔的乐高积木。

Joe和我一直在努力减少踩踏乐高积木的机会。一段时间后，我建议我们可以使用概率和统计数据来估计我们没捡到被孩子们落下的乐高积木的可能性。Joe非常赞同，“我的脚再也受不了踩到乐高积木上了！”。

我启动了我最喜欢的估算概率的工具，Joe和我开始研究在我们的扫描之后剩下乐高积木的可能性。

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')

np.random.seed(42)  # It's nice to replicate even virtual experiments

import pymc3 as pm
import scipy.stats as stats
print('Running on PyMC3 v{}'.format(pm.__version__))

> Running on PyMC3 v3.6

实验主义者做了二十次实验

 

每个人似乎都采用不同的方法来挑选乐高积木。对于Joe来说，他说他会在他的孩子们之后扫荡，然后拿起剩余的乐高积木。

我对Joe的第一个建议是，如果我们知道Joe在每次扫描中拾取乐高积木有多细致，那么我们可以确定每次扫描后留下乐高积木的可能性。数学上这是

 

其中p_Pickup Lego per sweep是在单次扫描中拾取乐高的概率，n_sweeps是Joe执行的扫描次数，而p_Lego remaining是完成所有扫描后Lego剩余的累积概率。

我建议Joe做一个实验，以确定他在拾取乐高积木方面有多细密：我会去他的孩子们玩乐高积木的房间，并在地板上散布随机数量的乐高积木，代表被孩子们随机扔下的乐高积木。乔将跟随我，捡起他发现的乐高积木。

我们需要重复这个实验多次，我估计20次，以便收敛到Joe对拾取乐高积木的有效性的估计。在这一点上，Joe说，“我想找一个解决方案，但是做二十次重复拾取乐高积木并不是我的乐趣。你需要向我证明这是值得我的时间的”。我同意这个实验似乎会伤害这个主题，所以我决定在我们进行实际实验之前向Joe证明这可能有用。

我设计了（虚拟）实验，为每个实验传播随机数量的乐高积木。基于Joe的一些估计，为了最好地代表Joe的孩子留下的乐高积木数量的实际情况，我选择了Legos的普通（又称高斯）分布，以平均100个乐高积木和20个乐高积木的标准偏差为中心对于每个实验：

 

其中N _0是Joe留下的Legos数量。

mean = 100
standard_deviation = 20
experiments = 20
N_0 = (np.random.randn(experiments)*standard_deviation+mean).astype(int)
N_0

> array([109,  97, 112, 130,  95,  95, 131, 115,  90, 110,  90,  90, 104, 61,  65,  88,  79, 106,  81,  71])

拿起乐高积木

 

然后，如果我们正在进行实际的实验，每次我散落完乐高积木之后，Joes（多次实验，所以是一个群体）会跟在我后面并拿起乐高积木。为了在虚拟实验中模拟这一点，我使用二项分布模拟了Legos Joe获取的数量。考虑到这一点，对于地面上的每个乐高，Joe会根据Joes拾取乐高的未知概率来拾取或不拾取它。

 

其中X是Joe发现的乐高积木的数量，而p是Joe拾取乐高积木的未知概率。这里也假设Joe在每次拾取积木时p始终不变。

X = np.random.binomial(N_0, actual_prob)  
X

> array([ 87,  66,  80, 110,  76,  65,  95,  92,  61,  83,  57,  76,  77, 46,  49,  75,  54,  75,  61,  54])

直方图显示了在虚拟实验的每个试验中已获得的乐高积木百分比的分布。

plt.hist(X / N_0, range=(0,1), bins=20)
plt.xlabel("Percentage Legos Picked Up")
plt.ylabel("Number of Times Observed")
plt.show();

 

对Joe进行建模

 

现在，我们用二项分布建了一个 Joe在每次扫描时拾起了多少个乐高积木的模型; 我们知道在每次实验中孩子们会留下多少乐高积木N _0，我们知道每次实验扫描中被拾起的乐高积木数目 X。我们不知道的是Joe拾起乐高积木的可能性，因此我们需要对概率进行建模。

概率分布的常见模型是β分布。在此示例中使用β分布有很多原因。β分布是二项分布的共轭先验; 这个原因是代数方面的，对于这个例子不太重要，因为我们将使用数字积分。更重要的是，β分布是用于将百分比或比例建模为随机变量的适当分布。

在我们的例子中，我们将使用一个弱信息的先验，一个估计Joe拾取乐高的概率在[0,1]的范围内，更接近该范围的中心。这说明我们对Joe在拾取乐高积木方面的技巧有所了解，并将其包含在模型中。β分布由两个值参数化，通常称为α（α）和β（β）。我选择的值与弱信息先验的目标相匹配，为0.5。

alpha_prior=2
beta_prior=2

x_locs = np.linspace(0, 1, 100)
plt.plot(x_locs, stats.beta.pdf(x_locs, alpha_prior, beta_prior), label='Probability Distribution Function');
plt.legend(loc='best');

 

模型拟合

现在我们有一个完整的模型可以生成靠观察才能获得的数据（即使到目前为止，这只是一个思想实验）：

 

使用PyMC3进行计算统计是很酷的。有很多文档和很好的介绍 ; 我不会在这里重复这些信息，而是直接跳到模型拟合。

首先我们构建模型对象：

basic_model = pm.Model()

with basic_model:
    p = pm.Beta('p', alpha=alpha_prior, beta=beta_prior)
    x = pm.Binomial('x', n=N_0, p=p, observed=X)

然后我将使用默认的No-U-Turn采样器（NUTS）来拟合模型：

with basic_model:
    trace_basic = pm.sample(50000, random_seed=123, progressbar=True)

> Auto-assigning NUTS sampler...
> Initializing NUTS using jitter+adapt_diag...
> Multiprocess sampling (2 chains in 2 jobs)
> NUTS: [p]
> Sampling 2 chains: 100%|██████████| 101000/101000 [00:55<00:00, 1828.39draws/s]

现在我们可以看一下模型拟合的一些结果。在这种情况下，我的朋友乔在一次捡拾中获得乐高的数据和给定模型（也称为后验概率）的估计为75％，95％的可能性置信区间为[73％， 77％]。我还绘制了先验概率（使用β分布的弱信息先验）与后验分布相比较的情况。

 

plt.hist(trace_basic['p'], 15, histtype='step', density=True, label='Posterior');
plt.plot(x_locs, stats.beta.pdf(x_locs, alpha_prior, beta_prior), label='Prior');
plt.legend(loc='best');

 

basic_model_df = pm.summary(trace_basic)
basic_model_df.round(2)

 

再次对乔建模

乔和我花了一些时间在这个模型上，看它是否适合我们对现实的感觉。毕竟，为了实现这一目标，乔需要进行大约二十次的实验，他希望有一些信心，这是值得他的时间。

假设Joe经常在95％置信区间的低端执行，我们学到的第一件事就是在经过4次扫描之后，剩下的任何乐高积木的可能性不到1％。 。

(1-basic_model_df.loc['p','hpd_2.5'])**4

> 0.0052992694812835335

整体而言，Joe和我对这个模型感到满意，但我们怀疑某些东西需要改进。Joe说他经常在拾取乐高积木时扫描房间四次，但是下次穿过房间时似乎总还会踩到乐高积木。我们更深入的探讨，我知道有时他会快速扫描，有时他会在房间里进行相当详细的扫描。

当我们最初模仿Joe时，我们认为Joe获得乐高的概率是一致的。所有统计模型都遗漏了一些东西，我们的原始模型也是如此。我们现在学到的是拾取一个离散的乐高的概率是不同的。现在，需要在模型中包含对生成数据的新理解。

有一个模型，称为Beta二项分布。Beta二项分布放宽了每个二项式试验存在单一概率的假设，建模每个二项式试验的概率参数不同。这与Joe描述的过程相匹配，有些扫描很快，有些非常详细。我们的新模型如下所示：

 

我们可以直接在PyMC3中对此进行建模。为此，我们提供半Cauchy分布作为β二项分布的α和β参数的弱正则化先验。我们使用半Cauchy分布作为一种方法来“鼓励”α和β值接近零，而不是将先验设置为均匀分布时发生的情况。在[0，∞）的范围内支持半Cauchy分布。有了它，我们有一个完整的

模型：

model_bb = pm.Model()

with model_bb:
    alpha_bb = pm.HalfCauchy('alpha_bb', beta = 2.)
    beta_bb = pm.HalfCauchy('beta_bb', beta = 2.)
    X_bb = pm.BetaBinomial('X_bb', alpha=alpha_bb, beta=beta_bb, n=N_0, observed=X)

with model_bb:
    trace_bb = pm.sample(50000, tuning=5000, random_seed=123, progressbar=True)

> Auto-assigning NUTS sampler...
> Initializing NUTS using jitter+adapt_diag...
> Multiprocess sampling (2 chains in 2 jobs)
> NUTS: [beta_bb, alpha_bb]
> Sampling 2 chains: 100%|██████████| 101000/101000 [03:05<00:00, 544.28draws/s]

通过这种新的参数化，我们失去了与概率参数的直接联系。这是必要的，因此Joe可以确定他需要完成多少次扫描才能达到他所需的水平的信心，即所有乐高积木都已被移除，并且他的脚在晚上可以安全地走在地板上。

在PyMC3中，我们可以通过基于拟合模型生成数据来确定总体概率后验估计。PyMC3使用后验预测检查有一种简单的方法。我将生成1000个后验概率的例子。

 

with model_bb:
    p_bb = pm.Beta('p_bb', alpha=alpha_bb, beta=beta_bb)
    ppc = pm.sample_posterior_predictive(trace_bb, 1000, vars=[p_bb])

> 100%|██████████| 1000/1000 [00:24<00:00, 41.64it/s]

有了这个，Joe和我使用β二项式假设（每个二项式试验的不同概率，或Legos的最低值的扫描）与二项式假设（所有二项式试验的单一概率）比较结果。当我们学习并预期β二项式模型假设中的概率分布比二项式假设更宽。

plt.hist(trace_basic['p'], 15, histtype='step', density=True, label='Posterior Binomial');
plt.hist(ppc['p_bb'], 15, histtype='step', density=True, label='Posterior BetaBinomial');
plt.plot(x_locs, stats.beta.pdf(x_locs, alpha_prior, beta_prior), label='Prior');
plt.legend(loc='best');

 

bb_quantiles = np.quantile(ppc['p_bb'], [0.025, 0.5, 0.975])
bb_quantiles

> array([0.59356599, 0.74900266, 0.86401046])

再次，做出安全的假设，Joe经常在95％置信区间的低端执行，我们首先得知的是，经过7次扫描后，剩下的任何乐高积木的可能性不到1％。已接晚上走路不会踩到积木的安全范围。

(1-bb_quantiles[0])**7

> 0.0018320202853132318

最后：模型与生成函数相比较

 

请记住，回到开头所有这一切都是一个思考练习，看看Joe是否值得进行20次实验，这样我们就可以确定Joe在扫描中获得乐高的概率。在这个思考练习中，我们生成的数据是，在20次实验扫描中，有多少乐高被搜集了。现在我们已经完成了建模，我们可以探索实际的数据生成功能并将其与模型进行比较。生成数据然后恢复参数，以验证建模方法的参数与原始参数有多接近的做法是计算统计中的最佳实践。

生成函数的这些参数在Jupyter笔记本中可用，位于开头附近的隐藏单元格中。

如下所示，在通过中拾取乐高的概率的原始生成函数是每次扫描乔所产生的β分布生成概率。由此可以看出，与二项式模型相比，beta二项式模型在生成数据中更好地重建原始生成函数。二项式模型没有考虑到Joe捡拾过程中的工作质量的变化，而beta二项模型确实如此。

plt.hist(actual_prob, label='Observed Probabilities')
plt.plot(x_locs, stats.beta.pdf(x_locs, alpha, beta), label='Generating Distribution');
plt.hist(trace_basic['p'], 15, histtype='step', density=True, label='Posterior Binomial');
plt.hist(ppc['p_bb'], 15, histtype='step', density=True, label='Posterior BetaBinomial');
plt.legend(loc='best');

 

乔不想做二十次

Joe和我对虚拟实验感到满意，并且相信通过执行实验，我们可以了解Joe在捡拾过程中获得乐高的概率。但Joe仍然不想做实验，“不是不想做二十次，而是一次都不想。我讨厌捡乐高积木，为什么我会一遍又一遍地学习一个关于我自己的参数。当然，汉克，必须有一个更好的方法。“

乔和我开始探索不同的方式，我们可以帮助他获得拾起所有乐高积木的信心，而无需对Joe进行实验。请继续关注下次。

这篇文章的完整Jupyter笔记本可用。