原文：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=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');

?

模型擬合

現在我們有一個完整的模型可以生成靠觀察才能獲得的數據（即使到目前為止，這只是一個思想實驗）：

?

使用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筆記本可用。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。