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

?

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

我有一個(gè)很好的朋友Joe，本周路過他家時(shí)我順便造訪了他家。像平常一樣，我們聊了天氣（在太平洋西北地區(qū)已經(jīng)比正常情況更熱），新聞（主要是關(guān)于我們?cè)撊绾螒?yīng)對(duì)炎熱天氣）和我們的孩子。我們倆都有孩子，而且他們都非常喜歡玩樂高積木。家里散落著大量的樂高積木就不可避免地會(huì)出現(xiàn)踩踏樂高積木的強(qiáng)烈痛苦，通常是在半夜或早上醒來制作咖啡時(shí)。為了避免出現(xiàn)那樣的痛苦，我們就追在孩子的后面收拾被他們亂扔的樂高積木。

Joe和我一直在努力減少踩踏樂高積木的機(jī)會(huì)。一段時(shí)間后，我建議我們可以使用概率和統(tǒng)計(jì)數(shù)據(jù)來估計(jì)我們沒撿到被孩子們落下的樂高積木的可能性。Joe非常贊同，“我的腳再也受不了踩到樂高積木上了！”。

我啟動(dòng)了我最喜歡的估算概率的工具，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

實(shí)驗(yàn)主義者做了二十次實(shí)驗(yàn)

?

每個(gè)人似乎都采用不同的方法來挑選樂高積木。對(duì)于Joe來說，他說他會(huì)在他的孩子們之后掃蕩，然后拿起剩余的樂高積木。

我對(duì)Joe的第一個(gè)建議是，如果我們知道Joe在每次掃描中拾取樂高積木有多細(xì)致，那么我們可以確定每次掃描后留下樂高積木的可能性。數(shù)學(xué)上這是

?

其中p_Pickup Lego per sweep是在單次掃描中拾取樂高的概率，n_sweeps是Joe執(zhí)行的掃描次數(shù)，而p_Lego remaining是完成所有掃描后Lego剩余的累積概率。

我建議Joe做一個(gè)實(shí)驗(yàn)，以確定他在拾取樂高積木方面有多細(xì)密：我會(huì)去他的孩子們玩樂高積木的房間，并在地板上散布隨機(jī)數(shù)量的樂高積木，代表被孩子們隨機(jī)扔下的樂高積木。喬將跟隨我，撿起他發(fā)現(xiàn)的樂高積木。

我們需要重復(fù)這個(gè)實(shí)驗(yàn)多次，我估計(jì)20次，以便收斂到Joe對(duì)拾取樂高積木的有效性的估計(jì)。在這一點(diǎn)上，Joe說，“我想找一個(gè)解決方案，但是做二十次重復(fù)拾取樂高積木并不是我的樂趣。你需要向我證明這是值得我的時(shí)間的”。我同意這個(gè)實(shí)驗(yàn)似乎會(huì)傷害這個(gè)主題，所以我決定在我們進(jìn)行實(shí)際實(shí)驗(yàn)之前向Joe證明這可能有用。

我設(shè)計(jì)了（虛擬）實(shí)驗(yàn)，為每個(gè)實(shí)驗(yàn)傳播隨機(jī)數(shù)量的樂高積木。基于Joe的一些估計(jì)，為了最好地代表Joe的孩子留下的樂高積木數(shù)量的實(shí)際情況，我選擇了Legos的普通（又稱高斯）分布，以平均100個(gè)樂高積木和20個(gè)樂高積木的標(biāo)準(zhǔn)偏差為中心對(duì)于每個(gè)實(shí)驗(yàn)：

?

其中N _0是Joe留下的Legos數(shù)量。

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])

拿起樂高積木

?

然后，如果我們正在進(jìn)行實(shí)際的實(shí)驗(yàn)，每次我散落完樂高積木之后，Joes（多次實(shí)驗(yàn)，所以是一個(gè)群體）會(huì)跟在我后面并拿起樂高積木。為了在虛擬實(shí)驗(yàn)中模擬這一點(diǎn)，我使用二項(xiàng)分布模擬了Legos Joe獲取的數(shù)量。考慮到這一點(diǎn)，對(duì)于地面上的每個(gè)樂高，Joe會(huì)根據(jù)Joes拾取樂高的未知概率來拾取或不拾取它。

?

其中X是Joe發(fā)現(xiàn)的樂高積木的數(shù)量，而p是Joe拾取樂高積木的未知概率。這里也假設(shè)Joe在每次拾取積木時(shí)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])

直方圖顯示了在虛擬實(shí)驗(yàn)的每個(gè)試驗(yàn)中已獲得的樂高積木百分比的分布。

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

?

對(duì)Joe進(jìn)行建模

?

現(xiàn)在，我們用二項(xiàng)分布建了一個(gè)?Joe在每次掃描時(shí)拾起了多少個(gè)樂高積木的模型;?我們知道在每次實(shí)驗(yàn)中孩子們會(huì)留下多少樂高積木N _0，我們知道每次實(shí)驗(yàn)掃描中被拾起的樂高積木數(shù)目?X。我們不知道的是Joe拾起樂高積木的可能性，因此我們需要對(duì)概率進(jìn)行建模。

概率分布的常見模型是β分布。在此示例中使用β分布有很多原因。β分布是二項(xiàng)分布的共軛先驗(yàn);?這個(gè)原因是代數(shù)方面的，對(duì)于這個(gè)例子不太重要，因?yàn)槲覀儗⑹褂脭?shù)字積分。更重要的是，β分布是用于將百分比或比例建模為隨機(jī)變量的適當(dāng)分布。

在我們的例子中，我們將使用一個(gè)弱信息的先驗(yàn)，一個(gè)估計(jì)Joe拾取樂高的概率在[0,1]的范圍內(nèi)，更接近該范圍的中心。這說明我們對(duì)Joe在拾取樂高積木方面的技巧有所了解，并將其包含在模型中。β分布由兩個(gè)值參數(shù)化，通常稱為α（α）和β（β）。我選擇的值與弱信息先驗(yàn)的目標(biāo)相匹配，為0.5。

alpha_prior=2 beta_prior=2x_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');

?

模型擬合

現(xiàn)在我們有一個(gè)完整的模型可以生成靠觀察才能獲得的數(shù)據(jù)（即使到目前為止，這只是一個(gè)思想實(shí)驗(yàn)）：

?

使用PyMC3進(jìn)行計(jì)算統(tǒng)計(jì)是很酷的。有很多文檔和很好的介紹?;?我不會(huì)在這里重復(fù)這些信息，而是直接跳到模型擬合。

首先我們構(gòu)建模型對(duì)象：

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)

然后我將使用默認(rèn)的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]

現(xiàn)在我們可以看一下模型擬合的一些結(jié)果。在這種情況下，我的朋友喬在一次撿拾中獲得樂高的數(shù)據(jù)和給定模型（也稱為后驗(yàn)概率）的估計(jì)為75％，95％的可能性置信區(qū)間為[73％， 77％]。我還繪制了先驗(yàn)概率（使用β分布的弱信息先驗(yàn)）與后驗(yàn)分布相比較的情況。

?

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)

?

再次對(duì)喬建模

喬和我花了一些時(shí)間在這個(gè)模型上，看它是否適合我們對(duì)現(xiàn)實(shí)的感覺。畢竟，為了實(shí)現(xiàn)這一目標(biāo)，喬需要進(jìn)行大約二十次的實(shí)驗(yàn)，他希望有一些信心，這是值得他的時(shí)間。

假設(shè)Joe經(jīng)常在95％置信區(qū)間的低端執(zhí)行，我們學(xué)到的第一件事就是在經(jīng)過4次掃描之后，剩下的任何樂高積木的可能性不到1％。 。

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

整體而言，Joe和我對(duì)這個(gè)模型感到滿意，但我們懷疑某些東西需要改進(jìn)。Joe說他經(jīng)常在拾取樂高積木時(shí)掃描房間四次，但是下次穿過房間時(shí)似乎總還會(huì)踩到樂高積木。我們更深入的探討，我知道有時(shí)他會(huì)快速掃描，有時(shí)他會(huì)在房間里進(jìn)行相當(dāng)詳細(xì)的掃描。

當(dāng)我們最初模仿Joe時(shí)，我們認(rèn)為Joe獲得樂高的概率是一致的。所有統(tǒng)計(jì)模型都遺漏了一些東西，我們的原始模型也是如此。我們現(xiàn)在學(xué)到的是拾取一個(gè)離散的樂高的概率是不同的。現(xiàn)在，需要在模型中包含對(duì)生成數(shù)據(jù)的新理解。

有一個(gè)模型，稱為Beta二項(xiàng)分布。Beta二項(xiàng)分布放寬了每個(gè)二項(xiàng)式試驗(yàn)存在單一概率的假設(shè)，建模每個(gè)二項(xiàng)式試驗(yàn)的概率參數(shù)不同。這與Joe描述的過程相匹配，有些掃描很快，有些非常詳細(xì)。我們的新模型如下所示：

?

我們可以直接在PyMC3中對(duì)此進(jìn)行建模。為此，我們提供半Cauchy分布作為β二項(xiàng)分布的α和β參數(shù)的弱正則化先驗(yàn)。我們使用半Cauchy分布作為一種方法來“鼓勵(lì)”α和β值接近零，而不是將先驗(yàn)設(shè)置為均勻分布時(shí)發(fā)生的情況。在[0，∞）的范圍內(nèi)支持半Cauchy分布。有了它，我們有一個(gè)完整的

模型：

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]

通過這種新的參數(shù)化，我們失去了與概率參數(shù)的直接聯(lián)系。這是必要的，因此Joe可以確定他需要完成多少次掃描才能達(dá)到他所需的水平的信心，即所有樂高積木都已被移除，并且他的腳在晚上可以安全地走在地板上。

在PyMC3中，我們可以通過基于擬合模型生成數(shù)據(jù)來確定總體概率后驗(yàn)估計(jì)。PyMC3使用后驗(yàn)預(yù)測(cè)檢查有一種簡(jiǎn)單的方法。我將生成1000個(gè)后驗(yàn)概率的例子。

?

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]

有了這個(gè)，Joe和我使用β二項(xiàng)式假設(shè)（每個(gè)二項(xiàng)式試驗(yàn)的不同概率，或Legos的最低值的掃描）與二項(xiàng)式假設(shè)（所有二項(xiàng)式試驗(yàn)的單一概率）比較結(jié)果。當(dāng)我們學(xué)習(xí)并預(yù)期β二項(xiàng)式模型假設(shè)中的概率分布比二項(xiàng)式假設(shè)更寬。

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])

再次，做出安全的假設(shè)，Joe經(jīng)常在95％置信區(qū)間的低端執(zhí)行，我們首先得知的是，經(jīng)過7次掃描后，剩下的任何樂高積木的可能性不到1％。已接晚上走路不會(huì)踩到積木的安全范圍。

(1-bb_quantiles[0])**7 >?0.0018320202853132318

最后：模型與生成函數(shù)相比較

?

請(qǐng)記住，回到開頭所有這一切都是一個(gè)思考練習(xí)，看看Joe是否值得進(jìn)行20次實(shí)驗(yàn)，這樣我們就可以確定Joe在掃描中獲得樂高的概率。在這個(gè)思考練習(xí)中，我們生成的數(shù)據(jù)是，在20次實(shí)驗(yàn)掃描中，有多少樂高被搜集了。現(xiàn)在我們已經(jīng)完成了建模，我們可以探索實(shí)際的數(shù)據(jù)生成功能并將其與模型進(jìn)行比較。生成數(shù)據(jù)然后恢復(fù)參數(shù)，以驗(yàn)證建模方法的參數(shù)與原始參數(shù)有多接近的做法是計(jì)算統(tǒng)計(jì)中的最佳實(shí)踐。

生成函數(shù)的這些參數(shù)在Jupyter筆記本中可用，位于開頭附近的隱藏單元格中。

如下所示，在通過中拾取樂高的概率的原始生成函數(shù)是每次掃描喬所產(chǎn)生的β分布生成概率。由此可以看出，與二項(xiàng)式模型相比，beta二項(xiàng)式模型在生成數(shù)據(jù)中更好地重建原始生成函數(shù)。二項(xiàng)式模型沒有考慮到Joe撿拾過程中的工作質(zhì)量的變化，而beta二項(xiàng)模型確實(shí)如此。

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和我對(duì)虛擬實(shí)驗(yàn)感到滿意，并且相信通過執(zhí)行實(shí)驗(yàn)，我們可以了解Joe在撿拾過程中獲得樂高的概率。但Joe仍然不想做實(shí)驗(yàn)，“不是不想做二十次，而是一次都不想。我討厭撿樂高積木，為什么我會(huì)一遍又一遍地學(xué)習(xí)一個(gè)關(guān)于我自己的參數(shù)。當(dāng)然，漢克，必須有一個(gè)更好的方法。“

喬和我開始探索不同的方式，我們可以幫助他獲得拾起所有樂高積木的信心，而無需對(duì)Joe進(jìn)行實(shí)驗(yàn)。請(qǐng)繼續(xù)關(guān)注下次。

這篇文章的完整Jupyter筆記本可用。

總結(jié)

以上是生活随笔為你收集整理的发现贝叶斯的乐高积木的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。