基于PyMC的貝葉斯建模實(shí)戰(zhàn)

即使對于一個(gè)非數(shù)據(jù)科學(xué)家來說，貝葉斯統(tǒng)計(jì)這個(gè)術(shù)語也已經(jīng)很流行了。你可能在大學(xué)期間把它作為必修課之一來學(xué)習(xí)，而沒有意識到貝葉斯統(tǒng)計(jì)有多么重要。事實(shí)上，貝葉斯統(tǒng)計(jì)不僅僅是一種特定的方法，甚至是一類方法；它是一種完全不同的統(tǒng)計(jì)分析范式。

為什么貝葉斯統(tǒng)計(jì)如此重要

貝葉斯統(tǒng)計(jì)為你提供了在新數(shù)據(jù)的證據(jù)中更新你的評估工具，這是一個(gè)在許多現(xiàn)實(shí)世界場景中常見的概念，如跟蹤大流行病，預(yù)測經(jīng)濟(jì)趨勢，或預(yù)測氣候變化。貝葉斯統(tǒng)計(jì)是許多較著名的統(tǒng)計(jì)模型的支柱，如高斯過程。

重要的是，學(xué)習(xí)貝葉斯統(tǒng)計(jì)原理可以成為你作為一個(gè)數(shù)據(jù)科學(xué)家的寶貴財(cái)富，因?yàn)樗o你一個(gè)全新的視角來解決具有真實(shí)世界動(dòng)態(tài)數(shù)據(jù)來源的新問題。

這篇文章將介紹貝葉斯統(tǒng)計(jì)的基本理論，以及如何在Python中實(shí)現(xiàn)一個(gè)簡單的貝葉斯模型。

目錄表：

01 什么是貝葉斯統(tǒng)計(jì)？

02 貝葉斯編程簡介

03 貝葉斯的工作流程

04 建立一個(gè)簡單的貝葉斯模型

閑話少說，進(jìn)入主題！讓我們開始介紹貝葉斯統(tǒng)計(jì)編程。

01 什么是貝葉斯統(tǒng)計(jì)？

你可能會(huì)在互聯(lián)網(wǎng)上的某個(gè)地方或在你的課堂上看到這個(gè)方程式。

如果你沒有，也不要擔(dān)心，因?yàn)槲覍⑾蚰愫喴榻B貝葉斯的基本原則以及該公式的工作原理。

關(guān)鍵術(shù)語

上述貝葉斯公式的組成部分一般被稱為概率聲明。例如，在下面的后驗(yàn)概率聲明中，該術(shù)語的意思是 "給定觀測值y，theta（θ）的概率是多少 "。

Theta（θ）是這里的未知數(shù)，被稱為我們所關(guān)心的參數(shù)。參數(shù)的不確定性遵循一個(gè)特定的概率分布，可以使用與數(shù)據(jù)相關(guān)的模型組合來估計(jì)有關(guān)參數(shù)。

上述貝葉斯統(tǒng)計(jì)表述也被稱為反概率，因?yàn)樗菑挠^察到參數(shù)開始的。換言之，貝葉斯統(tǒng)計(jì)試圖從數(shù)據(jù)（效果）中推斷出假設(shè)（原因），而不是用數(shù)據(jù)來接受/拒絕工作假設(shè)。

貝葉斯公式

那么，貝葉斯公式告訴我們什么呢？

后驗(yàn)概率是我們想知道的主要部分，因?yàn)門heta（θ）是我們感興趣的參數(shù)。

觀察的可能性僅僅意味著，在Theta（θ）的特定值下，數(shù)據(jù)y在現(xiàn)實(shí)世界中出現(xiàn)的可能性有多大。

先驗(yàn)概率是我們對Theta (θ)應(yīng)該是什么樣子的最佳猜測（例如，也許它遵循正態(tài)或高斯分布）。

歸一化常數(shù)只是一個(gè)系數(shù)常數(shù)，使整個(gè)方程積分為1（因?yàn)楦怕什荒艿陀?和高于1）。

現(xiàn)在我們已經(jīng)涵蓋了貝葉斯統(tǒng)計(jì)的基本理論，讓我們開始為即將到來的貝葉斯編程教程進(jìn)行設(shè)置。

02 貝葉斯編程介紹

安裝

首先，安裝PyMC3作為我們執(zhí)行貝葉斯統(tǒng)計(jì)編程的首選庫。

1. 推薦使用conda

conda install -c conda-forge pymc3

2. 也可使用pip

pip install pymc3

獲取數(shù)據(jù)

我們將使用描述美國家庭中氡氣(Radon)濃度的氡氣數(shù)據(jù)集。氡氣已被證明是非吸煙者患肺癌的最高預(yù)測因素之一，其濃度通常與房屋的整體條件（例如，是否有地下室，等等）有關(guān)。

首先，在你的筆記本或終端運(yùn)行以下命令：

!wget "https://raw.githubusercontent.com/fonnesbeck/mcmc_pydata_london_2019/master/data/radon.csv"

確保你的數(shù)據(jù)位于你的筆記本的同一目錄內(nèi)。

數(shù)據(jù)探索

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as npradon = pd.read_csv('./radon.csv', index_col=0)radon.head()

我們注意到，有29列描述了一個(gè)家庭中氡（Radon）的濃度。

讓我們畫一張圖，顯示 "ANOKA "的氡的對數(shù)濃度分布，用一條垂直線來說明對數(shù)濃度為1.1。

anoka_radon = radon.query('county=="ANOKA"').log_radon
sns.distplot(anoka_radon, bins=16)plt.axvline(1.1)

ANOKA地區(qū)氡氣對數(shù)濃度超過1.1的家庭比例似乎相當(dāng)大，這是一個(gè)令人擔(dān)憂的趨勢......

03 貝葉斯工作流

現(xiàn)在我們有了數(shù)據(jù)，讓我們進(jìn)行貝葉斯推斷。一般來說，這個(gè)過程可以分解為以下三個(gè)步驟。

第1步：指定一個(gè)概率模型

這是作為建模者要多做選擇的地方。你將需要為一切指定最可能的概率分布函數(shù)（例如，正態(tài)或高斯、考奇、二項(xiàng)式、t分布、F分布，等等）。

我所說的一切，是指包括未知參數(shù)、數(shù)據(jù)、協(xié)變量、缺失數(shù)據(jù)、預(yù)測在內(nèi)的一切。所以，用不同的分布函數(shù)做實(shí)驗(yàn)，看看在現(xiàn)實(shí)世界的場景中如何起效。

第2步：計(jì)算后驗(yàn)分布

現(xiàn)在你將計(jì)算這個(gè)概率項(xiàng)，給定貝葉斯方程右邊的所有項(xiàng)。

第3步：檢查你的模型

與其他ML模型一樣，評估你的模型是關(guān)鍵。回到第一步，檢查你的假設(shè)是否有意義。如果沒有，改變概率分布函數(shù)，并反復(fù)重申。

04 建立一個(gè)簡單的貝葉斯模型

現(xiàn)在，我將向你介紹一個(gè)簡單的編程練習(xí)來建立你的第一個(gè)貝葉斯模型。

第1步：定義一個(gè)貝葉斯模型

首先，讓我們定義我們的氡氣——貝葉斯模型，有兩個(gè)參數(shù)，平均值（μ-"miu"）和其偏差（σ-"sigma"）。這些參數(shù)（μ和σ）還需要通過選擇對應(yīng)的分布函數(shù)來建立模型（記住：我們必須為所有參數(shù)定義概率分布）。

對于這些，我們選擇的函數(shù)是正態(tài)/高斯分布（μ=0，σ=10）和均勻分布。你可以在模型的驗(yàn)證檢查中重新校準(zhǔn)這些值，如上面步驟3所述。

from pymc3 import Model, Normal, Uniformwith Model() as radon_model: 
  μ = Normal(’μ’, mu=0, sd=10) 
  σ = Uniform(’σ’, 0, 10)

下一步是用另一個(gè)概率分布來編譯radon_model本身。

**with** radon_model:
  dist = Normal('dist', mu=μ, sd=σ, observed=anoka_radon)

第2步：用數(shù)據(jù)進(jìn)行模型擬合

現(xiàn)在，我們需要用數(shù)據(jù)來擬合這個(gè)模型（即訓(xùn)練）。

from pymc3 import sample 
 
**with** radon_model: 
  samples = sample(1000, tune=1000, cores=2, random_seed=12)

讓我們畫出我們的參數(shù)μ在訓(xùn)練后的分布情況，同時(shí)畫出95%的置信線。

from arviz import plot_posterior

plot_posterior(samples, var_names=['μ'], ref_val=1.1)

好吧，看來1.1的對數(shù)濃度可能不是那么糟糕，因?yàn)樗窃诜植嫉奈捕耍ㄖ挥?.2%的樣品的對數(shù)濃度大于1.1）。