基於pgmpy的貝葉斯推理流程大白話解釋
pgmpy是github上的一個開源項目，在網站上他的簡介只有很簡單的一句話——pgmpy is a python library for working with Probabilistic Graphical Models.
用咱大俗話，它就是一個python庫可以推理“概率圖”模型，那什么是概率圖呢？很簡單請看下圖：
cancer model
這是一個很經典的圖，我們叫他用於診斷癌症的“貝葉斯模型”，別管他叫啥了，回到我們原來的問題——什么是概率圖，上面這個圖還不能稱為概率圖為什么呢？很明顯，缺少概率啊。那有同學說了，“你的概率跑哪去了？”。
別急啊，baby! 加上就是了，問題是該怎么加呢，我們可以以常識去思考，上面這個圖呢，中間這個圓圈是癌症，周圍四個圓圈分別是空氣污染、吸煙、x射線、呼吸障礙。再看圓圈之間的箭頭，我們是不是可以有個簡單的因果關系：某人得癌症的幾率，是不是受空氣污染、吸煙的影響，同時得了癌症的人，也會導致其受到x射線輻射治療，以及呈現呼吸障礙的症狀啊。
這樣通過以上的因果分析，概率就明了啦，我們就可以去添加概率了，首先每個圓圈，我們再學術一點就稱為“節點”吧，對於每個節點，我們不考慮他們之間的關聯，僅在一個人群里面去觀察，人群中有多少人到了癌症，多少人沒得，這是不是有個得癌症的概率呢，學術一點，我們稱為“先驗概率”，同理其他的每個節點的先驗概率，我們都可以通過這種思想得到。
如果我們繼續觀察患癌人群，那么該人群中有呼吸障礙症狀的人的概率，我們是不是就知道了，得有同學說了——這也是“先驗概率”嗎？哈哈，當然不是，就看我的表述，是不是比前面的那個得到癌症概率的表述多了一個條件呢，我這里可是加了一個前提——患癌人群中有呼吸症狀的人的概率，所以我們稱之為"條件概率"。
我們再回到開始"概率圖"這個問題，這個概率到底指的是"先驗概率"還是“條件概率”呢，不科學的講，指的是條件概率，具體為啥說“不科學呢”，那是因為對於最上面的節點來說，沒有別的節點指向他們，那這個時候他們的先驗概率也是概率圖中的概率了。
到現在為止我們介紹完概率圖了，核心來了，啥叫人工智能啊。這就要對應我們文章開頭說到的推理了，所謂的推理，在這里我們稱為貝葉斯推理，有了前面的貝葉斯模型，再進行貝葉斯推理當然是水到渠成的事了，貝葉斯推理的目標是是什么呢，目標還是"概率"。說到這，可能有同學懵逼了，咋還是概率啊？哈哈，這也沒辦法，貝葉斯推理的核心也就是任何事情都是不確定性的表達，這中間還夾雜着傳統的頻率學派與貝葉斯學派的爭論，在這里我們不做深入的探討，有興趣的同學可以自行百度。
回到我們的正題，前面說到我們貝葉斯推理的目的還是個“概率”，那這個概率是什么呢？先說稱呼，學術上稱為“后驗概率”。舉個例子，我們有大量的數據表明：一個人患癌症與吸煙呈現明顯的關系，甚至我們已經得到兩者之間的條件概率，這時候呢，來了一個患者，他有吸煙的習慣還有點輕微的呼吸障礙的症狀，那我們想確定他患癌症的概率是多少？這個概率就是我們所求的“后驗概率”，而他抽煙的習慣以及呼吸障礙的症狀，都是我們事先掌握的”證據信息”，這個證據信息將在推理過程中起到關鍵作用，也是能得到我們該病人患癌症概率的關鍵。
下面我們用代碼展示了從建模到添加概率再到最終推理的實現過程，短短幾行代碼，就讓你學會應用人工智能的貝葉斯推理技術，是不是很炫酷呢~~
## step1 定義模型結構 from pgmpy.models import BayesianModel cancer_model = BayesianModel([('Pollution', 'Cancer'), ('Smoker', 'Cancer'), ('Cancer', 'Xray'), ('Cancer', 'Dyspnoea')]) ##step2 添加概率 from pgmpy.factors.discrete import TabularCPD cpd_poll = TabularCPD(variable='Pollution', variable_card=2, values=[[0.9], [0.1]]) cpd_smoke = TabularCPD(variable='Smoker', variable_card=2, values=[[0.3], [0.7]]) cpd_cancer = TabularCPD(variable='Cancer', variable_card=2, values=[[0.03, 0.05, 0.001, 0.02], [0.97, 0.95, 0.999, 0.98]], evidence=['Smoker', 'Pollution'], evidence_card=[2, 2]) cpd_xray = TabularCPD(variable='Xray', variable_card=2, values=[[0.9, 0.2], [0.1, 0.8]], evidence=['Cancer'], evidence_card=[2]) cpd_dysp = TabularCPD(variable='Dyspnoea', variable_card=2, values=[[0.65, 0.3], [0.35, 0.7]], evidence=['Cancer'], evidence_card=[2]) cancer_model.add_cpds(cpd_poll, cpd_smoke, cpd_cancer, cpd_xray, cpd_dysp) ##step3 進行推理 from pgmpy.inference import VariableElimination cancer_infer = VariableElimination(cancer_model) q = asia_infer.query(variables=['cancer'], evidence={'smoke': 'yes'}) print(q)pgmpy 教程
Table of Contents
創建貝葉斯網
在 pgmpy 中, 定義一個貝葉斯網的流程一般是先建立網絡結構, 然后填入相關參數.
建立網絡結構
from pgmpy.models import BayesianModel
cancer_model = BayesianModel([('Pollution', 'Cancer'),
('Smoker', 'Cancer'),
('Cancer', 'Xray'),
('Cancer', 'Dyspnoea')])
這個網絡中有五個節點: Pollution, Cancer, Smoker, Xray, Dyspnoea.
('Pollution', 'Cancer'): 一條有向邊, 從 Pollution 指向 Cancer, 表示環境污染有可能導致癌症.
('Smoker', 'Cancer'): 吸煙有可能導致癌症.
('Cancer', 'Xray'): 得癌症的人可能會去照X射線.
('Cancer', 'Dyspnoea'): 得癌症的人可能會呼吸困難.
填入參數
from pgmpy.factors.discrete import TabularCPD
cpd_poll = TabularCPD(variable='Pollution', variable_card=2,
values=[[0.9], [0.1]])
cpd_smoke = TabularCPD(variable='Smoker', variable_card=2,
values=[[0.3], [0.7]])
cpd_cancer = TabularCPD(variable='Cancer', variable_card=2,
values=[[0.03, 0.05, 0.001, 0.02],
[0.97, 0.95, 0.999, 0.98]],
evidence=['Smoker', 'Pollution'],
evidence_card=[2, 2])
cpd_xray = TabularCPD(variable='Xray', variable_card=2,
values=[[0.9, 0.2], [0.1, 0.8]],
evidence=['Cancer'], evidence_card=[2])
cpd_dysp = TabularCPD(variable='Dyspnoea', variable_card=2,
values=[[0.65, 0.3], [0.35, 0.7]],
evidence=['Cancer'], evidence_card=[2])
這些代碼, 主要是建立一些概率表, 然后往表里面填入了一些參數.
Pollution: 有兩種概率, 分別是 0.9 和 0.1.
Smoker: 有兩種概率, 分別是 0.3 和 0.7. (意思是在一個人群里, 有 30% 的人吸煙, 有 70% 的人不吸煙)
Cancer: envidence 表示有 Smoker 和 Pollution 兩個節點指向 Cancer 節點;
[0.03, 0.05, 0.001, 0.02] 的意思, 通過下表可以較容易看出來:
| Pollution | Pollution_0 | Pollution_0 | Pollution_1 | Pollution_1 |
|---|---|---|---|---|
| Smoker | Smoker_0 | Smoker_1 | Smoker_0 | Smoker_1 |
| Cancer_0 | 0.03 | 0.05 | 0.001 | 0.02 |
| Cancer_1 | 0.97 | 0.95 | 0.999 | 0.98 |
可以看出, 環境不好, 又抽煙的人, 有 0.98 的概率得癌症.
Xray: 通過 envidence 可以看出, 只有 Cancer 指向它; [0.9, 0.2] 表示, 未得癌症的人去照X光, 有 0.9 的概率正常, 而得了癌症的人去照X光, 有 0.2 的概率正常, [0.1, 0.8] 表示, 未得癌症的人去照X光, 有 0.1 的概率不正常, 而得癌症的人有 0.8 的概率不正常.
Dyspnoea: 未得癌症的人有 0.65 的可能性呼吸正常, 得癌症的人只有 0.3 的可能性呼吸正常; 未得癌症的人有 0.35 的可能性呼吸困難, 得癌症的人有 0.7 的可能性呼吸困難.
將數據加入網絡
# Associating the parameters with the model structure. cancer_model.add_cpds(cpd_poll, cpd_smoke, cpd_cancer, cpd_xray, cpd_dysp) # Checking if the cpds are valid for the model. cancer_model.check_model()
第一行代碼表示將概率分布表加入網絡中; 第二行代碼驗證模型數據的正確性(至於怎么驗證的, 我沒有看源碼, 但是猜測是驗證相關數據的和是否為1, 驗證數據中是否有負數, 等等).
相關的檢查
cancer_model.is_active_trail('Pollution', 'Smoker') # False
cancer_model.is_active_trail('Pollution', 'Smoker', observed=['Cancer']) # True
cancer_model.local_independencies('Xray') # (Xray _|_ Dyspnoea, Smoker, Pollution | Cancer)
cancer_model.get_independencies()
第一行代碼表示, Pollution 節點是否指向 Smoker 節點, 返回的是 False.
第二行代碼表示, Cancer 節點的計算, 是否與 Pollution 和 Smoker 兩個節點有關, 返回的是 True.
第三行代碼的輸出表示, Xray 與 Dyspnoea, Smoker, Pollution 不相關, 與 Cancer 相關.
第四行代碼輸出所有的依賴關系.
使用 Asia 網絡進行推理
在構建了貝葉斯網之后, 我們使用貝葉斯網來進行推理. 推理算法分精確推理和近似推理. 精確推理有變量消元法和團樹傳播法; 近似推理算法是基於隨機抽樣的算法.
Asia 網絡是早期貝葉斯網文獻中給出的一個網絡, 與肺結核, 肺癌, 支氣管炎等有關, 和我們之前的那個網絡很相似. 我們使用 Asia 網絡來看一下如何進行推理.
創建源文件 asia.py
下載網絡模型放到源文件目錄下:
wget http://www.bnlearn.com/bnrepository/asia/asia.bif.gz gzip -qd asia.bif.gz
編輯 asia.py:
導入 Asia 網絡:
from pgmpy.readwrite import BIFReader
reader = BIFReader('asia.bif')
asia_model = reader.get_model()
輸出節點信息:
print(asia_model.nodes()) # ['xray', 'bronc', 'asia', 'dysp', 'lung', 'either', 'smoke', 'tub']
輸出依賴關系:
print(asia_model.edges())
結果如下:
[('bronc', 'dysp'),
('asia', 'tub'),
('lung', 'either'),
('either', 'dysp'),
('either', 'xray'),
('smoke', 'bronc'),
('smoke', 'lung'),
('tub', 'either')]
查看概率分布:
print(asia_model.get_cpds())
結果如下:
[<TabularCPD representing P(asia:2) at 0x7f92cf1e46d0>, <TabularCPD representing P(bronc:2 | smoke:2) at 0x7f92ea70f990>, <TabularCPD representing P(dysp:2 | bronc:2, either:2) at 0x7f92f80f4ed0>, <TabularCPD representing P(either:2 | lung:2, tub:2) at 0x7f92e8e71dd0>, <TabularCPD representing P(lung:2 | smoke:2) at 0x7f92cf1e4c50>, <TabularCPD representing P(smoke:2) at 0x7f92cf1e4ed0>, <TabularCPD representing P(tub:2 | asia:2) at 0x7f92cf1ec210>, <TabularCPD representing P(xray:2 | either:2) at 0x7f92cf1ec410>]
直觀顯示某節點的概率分布:
print(asia_model.get_cpds('xray').values)
結果如下:
[[0.98 0.05] [0.02 0.95]]
變量消除法
變量消除法是精確推斷的一種方法.
from pgmpy.inference import VariableElimination
asia_infer = VariableElimination(asia_model)
q = asia_infer.query(variables=['bronc'], evidence={'smoke': 0})
print(q['bronc'])
結果如下:
+---------+--------------+ | bronc | phi(bronc) | |---------+--------------| | bronc_0 | 0.6000 | | bronc_1 | 0.4000 | +---------+--------------+
意思是, 在不吸煙的情況下, 得支氣管炎的概率是 0.4, 未得支氣管炎的概率是 0.6.
在在吸煙情況下, 得支氣管炎的概率和未得支氣管炎的概率可以這樣查詢:
q = asia_infer.query(variables=['bronc'], evidence={'smoke': 1})
print(q['bronc'])
結果如下:
+---------+--------------+ | bronc | phi(bronc) | |---------+--------------| | bronc_0 | 0.3000 | | bronc_1 | 0.7000 | +---------+--------------+
利用訓練數據學習
場景: 投硬幣, 訓練數據中, 有 30% 的正面朝上, 有 70% 的反面朝上. 我們使用極大似然估計和狄利克雷分布下貝葉斯參數先驗估計硬幣的條件概率分布.
生成數據:
import numpy as np import pandas as pd raw_data = np.array([0] * 30 + [1] * 70) data = pd.DataFrame(raw_data, columns=['coin'])
定義貝葉斯網:
from pgmpy.models import BayesianModel
from pgmpy.estimators import MaximumLikelihoodEstimator, BayesianEstimator
model = BayesianModel()
model.add_node('coin')
使用極大似然估計:
model.fit(data, estimator=MaximumLikelihoodEstimator)
print(model.get_cpds('coin'))
結果如下:
+---------+-----+ | coin(0) | 0.3 | +---------+-----+ | coin(1) | 0.7 | +---------+-----+
使用狄利克雷分布作為先驗分布的貝葉斯推論:
model.fit(data, estimator=BayesianEstimator, prior_type='dirichlet', pseudo_counts={'coin': [50, 50]})
print(model.get_cpds('coin'))
結果如下:
+---------+-----+ | coin(0) | 0.4 | +---------+-----+ | coin(1) | 0.6 | +---------+-----+
更復雜的例子(學生例子)
學生例子中包含5個隨機變量, 如下所示:
| 變量 | 含義 | 取值 |
|---|---|---|
| Difficulty | 課程本身難度 | 0=easy, 1=difficult |
| Intelligence | 學生聰明程度 | 0=stupid, 1=smart |
| Grade | 學生課程成績 | 1=A, 2=B, 3=C |
| SAT | 學生高考成績 | 0=low, 1=high |
| Letter | 可否獲得推薦信 | 0=未獲得, 1=獲得 |
生成數據:
import numpy as np import pandas as pd raw_data = np.random.randint(low=0, high=2, size=(1000, 5)) data = pd.DataFrame(raw_data, columns=['D', 'I', 'G', 'L', 'S'])
定義模型:
from pgmpy.models import BayesianModel
from pgmpy.estimators import MaximumLikelihoodEstimator, BayesianEstimator
model = BayesianModel([('D', 'G'), ('I', 'G'), ('I', 'S'), ('G', 'L')])
model.fit(data, estimator=MaximumLikelihoodEstimator)
for cpd in model.get_cpds():
print("CPD of {variable}:".format(variable=cpd.variable))
print(cpd)
結果如下:
CPD of D: +------+-------+ | D(0) | 0.501 | +------+-------+ | D(1) | 0.499 | +------+-------+ CPD of G: +------+------+----------------+----------------+----------------+ | D | D(0) | D(0) | D(1) | D(1) | +------+------+----------------+----------------+----------------+ | I | I(0) | I(1) | I(0) | I(1) | +------+------+----------------+----------------+----------------+ | G(0) | 0.48 | 0.509960159363 | 0.444915254237 | 0.551330798479 | +------+------+----------------+----------------+----------------+ | G(1) | 0.52 | 0.490039840637 | 0.555084745763 | 0.448669201521 | +------+------+----------------+----------------+----------------+ CPD of I: +------+-------+ | I(0) | 0.486 | +------+-------+ | I(1) | 0.514 | +------+-------+ CPD of L: +------+----------------+----------------+ | G | G(0) | G(1) | +------+----------------+----------------+ | L(0) | 0.489959839357 | 0.501992031873 | +------+----------------+----------------+ | L(1) | 0.510040160643 | 0.498007968127 | +------+----------------+----------------+ CPD of S: +------+----------------+----------------+ | I | I(0) | I(1) | +------+----------------+----------------+ | S(0) | 0.512345679012 | 0.468871595331 | +------+----------------+----------------+ | S(1) | 0.487654320988 | 0.531128404669 | +------+----------------+----------------+
參數學習:
pseudo_counts = {'D': [300, 700], 'I': [500, 500], 'G': [800, 200], 'L': [500, 500], 'S': [400, 600]}
model.fit(data, estimator=BayesianEstimator, prior_type='dirichlet', pseudo_counts=pseudo_counts)
for cpd in model.get_cpds():
print("CPD of {variable}:".format(variable=cpd.variable))
print(cpd)
結果如下:
CPD of D: +------+--------+ | D(0) | 0.4005 | +------+--------+ | D(1) | 0.5995 | +------+--------+ CPD of G: +------+-------+----------------+----------------+----------------+ | D | D(0) | D(0) | D(1) | D(1) | +------+-------+----------------+----------------+----------------+ | I | I(0) | I(1) | I(0) | I(1) | +------+-------+----------------+----------------+----------------+ | G(0) | 0.736 | 0.741806554756 | 0.732200647249 | 0.748218527316 | +------+-------+----------------+----------------+----------------+ | G(1) | 0.264 | 0.258193445244 | 0.267799352751 | 0.251781472684 | +------+-------+----------------+----------------+----------------+ CPD of I: +------+-------+ | I(0) | 0.493 | +------+-------+ | I(1) | 0.507 | +------+-------+ CPD of L: +------+----------------+----------------+ | G | G(0) | G(1) | +------+----------------+----------------+ | L(0) | 0.496662216288 | 0.500665778961 | +------+----------------+----------------+ | L(1) | 0.503337783712 | 0.499334221039 | +------+----------------+----------------+ CPD of S: +------+----------------+----------------+ | I | I(0) | I(1) | +------+----------------+----------------+ | S(0) | 0.436742934051 | 0.423381770145 | +------+----------------+----------------+ | S(1) | 0.563257065949 | 0.576618229855 | +------+----------------+----------------+