貝葉斯網絡python實戰(以泰坦尼克號數據集為例,pgmpy庫)


貝葉斯網絡python實戰(以泰坦尼克號數據集為例,pgmpy庫)

 


本文的相關數據集,代碼見文末百度雲

 

貝葉斯網絡簡介

貝葉斯網絡是一種置信網絡,一個生成模型。(判別模型,生成模型的區分可以這樣:回答p(label|x)即樣本x屬於某一類別的可能的,就是判別模型,而回答p(x,label) 和p(x|label)的,即回答在給定的類別中找樣本x及樣本分布情況的,即為生成模型。生成模型給出的聯合分布比判別網絡能給出更多的信息,對其求邊緣分布即可得p(label|x) p(x|label))同時貝葉斯網絡還是一個簡單的白盒網絡,提供了高可解釋性的可能。相比於大熱的幾乎無所不能的深度神經網絡,貝葉斯網絡仍有他的優勢和應用場景。比如在故障分析,疾病診斷里,我們不僅需要回答是不是,更重要的是回答為什么,並給出依據。這樣的場景下,以貝葉斯網絡為代表的一些可解釋好的白盒網絡更加有優勢。

貝葉斯推斷思路

與頻率派直接從數據統計分析構建模型不同,貝葉斯派引入一個先驗概率,表達對事件的已有了解,然后利用觀測數據對先驗知識進行修正,如通常把拋硬幣向上的概率認為是0.5,這是個很朴素的先驗知識,若是實驗結果拋出了500-500的結果,那么證明先驗知識是可靠的,合適的,若是出現100-900結果,那么先驗知識會被逐漸修改(越來越相信這是個作弊硬幣),當實驗數據足夠多的時候,先驗知識就幾乎不再體現,這時候得到與頻率派幾乎相同的結果。如圖
在這里插入圖片描述

具體例子推導可見here

貝葉斯網絡

貝葉斯網絡結構如下所示,其是有特征節點和鏈接構成的有向無環圖。節點上是概率P(A),P(B)… 連接上是條件概率P(A|B) P(A|C) … 即若有A指向B的連接,則連接代表的就應為P(B|A),更多信息可參考以下內容,這里不再贅述,貝葉斯網絡結構本身不困難,其難點主要在於推理算法等數值計算問題,如為應用則無需深究。
貝葉斯網絡發展及其應用綜述
《貝葉斯網絡引論》@張連文
靜態貝葉斯網絡
在這里插入圖片描述

貝葉斯網絡的實現

相關工具一直很豐富,matlab,R上都有成熟的工具。這里使用了python下的pgmpy,輕量好用,不像pymc那樣容易安裝困難。
安裝:
conda install -c ankurankan pgmpy

pip install pgmpy

應用步驟

1.先確定以那些變量(特征)為節點,這里還包括由特征工程特征選擇之類的工作。當然若有專業知識的參與會得到更合理的特征選擇。
2.確定網絡結構(拓撲)用以反應變量節點之間的依賴關系。也就是明確圖的結構。這里既可以在有專家參與的情況下手工設計,也可以自動找到高效合適的網絡,稱為結構學習。貝葉斯網絡的結構對最終網絡性能很關鍵,若是構建所謂全連接貝葉斯網(即各個變量間兩兩相連),雖沒有遺漏關聯,但會導致嚴重的過擬合,因為數據量很難支撐起全連接直接海量的條件概率。
3.明確每條邊上的條件概率。和結構一樣,參數也可由專家手工確定(先驗),亦可通過數據自動學習(即參數學習),或兩者同時進行。
下面以一個經典數據集為例展示如何利用pgmpy包進行貝葉斯網絡建模

泰坦尼克數據集背景介紹

ref:https://www.jianshu.com/p/9b6ee1fb7a60
https://www.kaggle.com/c/titanic
這是kaggle經典數據集,主要是讓參賽選手根據訓練集中的乘客數據和存活情況進行建模,進而使用模型預測測試集中的乘客是否會存活。乘客特征總共有11個,以下列出。這個數據集特征明確,數據量不大,很適合應用貝葉斯網絡之類的模型來做,目前最好的結果是正確率應該有80+%(具體多少因為答案泄露不好講了)

PassengerId => 乘客ID
Pclass => 客艙等級(1/2/3等艙位)
Name => 乘客姓名
Sex => 性別
Age => 年齡
SibSp => 兄弟姐妹數/配偶數
Parch => 父母數/子女數
Ticket => 船票編號
Fare => 船票價格
Cabin => 客艙號
Embarked => 登船港口
在開始建模之前,先進行下特征工程,處理原始數據集的缺項等。這里前面處理主要采用https://www.jianshu.com/p/9b6ee1fb7a60的方法(他應用pandas清理數據的技巧很值得一學),我在他的處理后,進一步進行了一些離散化處理,以使得數據符合貝葉斯網絡的要求(貝葉斯網絡也有支持連續變量的版本,但因為推理,學習的困難,目前還用的很少),最后保留5個特征。

''' PassengerId => 乘客ID Pclass => 客艙等級(1/2/3等艙位) Name => 乘客姓名 Sex => 性別 清洗成male=1 female=0 Age => 年齡 插補后分0,1,2 代表 幼年(0-15) 成年(15-55) 老年(55-) SibSp => 兄弟姐妹數/配偶數 Parch => 父母數/子女數 Ticket => 船票編號 Fare => 船票價格 經聚類變0 1 2 代表少 多 很多 Cabin => 客艙號 清洗成有無此項,並發現有的生存率高 Embarked => 登船港口 清洗na,填S ''' # combine train and test set. train=pd.read_csv('./train.csv') test=pd.read_csv('./test.csv') full=pd.concat([train,test],ignore_index=True) full['Embarked'].fillna('S',inplace=True) full.Fare.fillna(full[full.Pclass==3]['Fare'].median(),inplace=True) full.loc[full.Cabin.notnull(),'Cabin']=1 full.loc[full.Cabin.isnull(),'Cabin']=0 full.loc[full['Sex']=='male','Sex']=1 full.loc[full['Sex']=='female','Sex']=0 full['Title']=full['Name'].apply(lambda x: x.split(',')[1].split('.')[0].strip()) nn={'Capt':'Rareman', 'Col':'Rareman','Don':'Rareman','Dona':'Rarewoman', 'Dr':'Rareman','Jonkheer':'Rareman','Lady':'Rarewoman','Major':'Rareman', 'Master':'Master','Miss':'Miss','Mlle':'Rarewoman','Mme':'Rarewoman', 'Mr':'Mr','Mrs':'Mrs','Ms':'Rarewoman','Rev':'Mr','Sir':'Rareman', 'the Countess':'Rarewoman'} full.Title=full.Title.map(nn) # assign the female 'Dr' to 'Rarewoman' full.loc[full.PassengerId==797,'Title']='Rarewoman' full.Age.fillna(999,inplace=True) def girl(aa): if (aa.Age!=999)&(aa.Title=='Miss')&(aa.Age<=14): return 'Girl' elif (aa.Age==999)&(aa.Title=='Miss')&(aa.Parch!=0): return 'Girl' else: return aa.Title full['Title']=full.apply(girl,axis=1) Tit=['Mr','Miss','Mrs','Master','Girl','Rareman','Rarewoman'] for i in Tit: full.loc[(full.Age==999)&(full.Title==i),'Age']=full.loc[full.Title==i,'Age'].median() full.loc[full['Age']<=15,'Age']=0 full.loc[(full['Age']>15)&(full['Age']<55),'Age']=1 full.loc[full['Age']>=55,'Age']=2 full['Pclass']=full['Pclass']-1 from sklearn.cluster import KMeans Fare=full['Fare'].values Fare=Fare.reshape(-1,1) km = KMeans(n_clusters=3).fit(Fare) #將數據集分為2類 Fare = km.fit_predict(Fare) full['Fare']=Fare full['Fare']=full['Fare'].astype(int) full['Age']=full['Age'].astype(int) full['Cabin']=full['Cabin'].astype(int) full['Pclass']=full['Pclass'].astype(int) full['Sex']=full['Sex'].astype(int) #full['Survived']=full['Survived'].astype(int) dataset=full.drop(columns=['Embarked','Name','Parch','PassengerId','SibSp','Ticket','Title']) dataset.dropna(inplace=True) dataset['Survived']=dataset['Survived'].astype(int) #dataset=pd.concat([dataset, pd.DataFrame(columns=['Pri'])]) train=dataset[:800] test=dataset[800:] ''' 最后保留如下項目,並切出800的訓練集: Pclass => 客艙等級(0/1/2等艙位) Sex => 性別 male=1 female=0 Age => 年齡 插補后分0,1,2 代表 幼年(0-15) 成年(15-55) 老年(55-) Fare => 船票價格 經聚類變0 1 2 代表少 多 很多 Cabin => 客艙號 清洗成有無此項,並發現有的生存率高 ''' 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80

模型結構搭建

這里先手動設計網絡結構。
憑借對數據的理解,先設計如下的結構在這里插入圖片描述

from pgmpy.models import BayesianModel from pgmpy.estimators import BayesianEstimator #model = BayesianModel([('Age', 'Pri'), ('Sex', 'Pri'),('Pri','Survived'),('Fare','Pclass'),('Pclass','Survived'),('Cabin','Survived')]) model = BayesianModel([('Age', 'Survived'), ('Sex', 'Survived'),('Fare','Pclass'),('Pclass','Survived'),('Cabin','Survived')]) 
  • 1
  • 2
  • 3
  • 4
  • 5

其中(‘Age’, ‘Survived’)表示Age指向Survived

pgmpy沒有提供可視化,這里簡單用graphviz實現了一下。

def showBN(model,save=False): '''傳入BayesianModel對象,調用graphviz繪制結構圖,jupyter中可直接顯示''' from graphviz import Digraph node_attr = dict( style='filled', shape='box', align='left', fontsize='12', ranksep='0.1', height='0.2' ) dot = Digraph(node_attr=node_attr, graph_attr=dict(size="12,12")) seen = set() edges=model.edges() for a,b in edges: dot.edge(a,b) if save: dot.view(cleanup=True) return dot showBN(model) 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20

模型參數構建

接下來就是確定網絡的參數,也就是各個邊上的條件概率。
若手工填入,可這樣寫
在這里插入圖片描述

from pgmpy.factors.discrete import TabularCPD #構建cpd表 my_cpd= TabularCPD(variable='Pclass', variable_card=3, values=[[0.65, 0.3], [0.30, 0.6],[0.05,0.1]], evidence=['Fare'], evidence_card=[2]) # 填cpd表 model.add_cpds(my_cpd) # 執行檢查(可選,用於檢查cpd是否填錯) cancer_model.check_model() 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

但在此例子里,使用參數學習的辦法從數據里自動學習。
參數學習有兩種典型方法,極大似然估計和貝葉斯估計。因為前者的過擬合嚴重,一般都使用后者進行參數學習。pgmpy提供的貝葉斯估計器提供三種先驗分布的支持,‘dirichlet’, ‘BDeu’, ‘K2’,實際上都是dirichlet分布,這里解釋下貝葉斯估計器的工作原理。

貝葉斯估計器

在貝葉斯分析的框架下,待求參數θ被看做是隨機變量,對他的估計就是在其先驗上,用數據求后驗,因此先要有對θ的先驗假設。而我們通常取的先驗分布就是dirichlet(狄利克雷)分布。對於一個含有i個離散狀態的節點,我們設其參數為在這里插入圖片描述
並令其先驗為狄利克雷分布D[α1,α2…αi](i=2時也稱beta分布)

在這里插入圖片描述
這個先驗有i個參數,數學上證明了,這些參數就相當於將先驗表達成了α個虛擬樣本,其中滿足X=xi的樣本數為αi,這個α就成為等價樣本量(equivalent_sample_size)。這個巧合其實正式先驗函數取這個函數的緣由,另外,其計算后的后驗分布也是狄利克雷分布(稱這種行為叫共軛先驗)。注:各個分布可參考 https://zhuanlan.zhihu.com/p/37976562
至此可以理解pgmpy提供的貝葉斯估計器的參數的含義了。
其定義為estimate_cpd(node, prior_type=‘BDeu’, pseudo_counts=[], equivalent_sample_size=5)
node是節點名
當prior_type=‘BDeu’ 就表示選擇了一個equivalent_sample_size=5的無差別客觀先驗,認定各個概率相等,不提供信息,但並不是沒用,這個先驗起到了類似神經網絡里頭控制過擬合的正則項的作用。
當prior_type='dirichlet’表示選擇一般的狄利克雷分布,這時候要主動填入[α1,α2…αi]
當prior_type= ‘K2’ 意為 ‘dirichlet’ + setting every pseudo_count to 1
具體使用如下:

data = pd.DataFrame(data={'A': [0, 0, 1], 'B': [0, 1, 0], 'C': [1, 1, 0]}) model = BayesianModel([('A', 'C'), ('B', 'C')]) estimator = BayesianEstimator(model, data) cpd_C = estimator.estimate_cpd('C', prior_type="dirichlet", pseudo_counts=[1, 2]) model.add_cpds(cpd_C) 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

上面是一個一個填進去的,在本例中有更簡單的方法,就是利用提供的fit函數,一並估計各個cpd(條件概率),即

model.fit(train, estimator=BayesianEstimator, prior_type="BDeu") # default equivalent_sample_size=5 
  • 1

直接把前面處理得到dataframe 傳入即可。這里記錄一個bug:pgmpy目前將離散變量命名限制為從0開始,所以本例子里的Pclass 項從(1/2/3等級)都減一處理成了(0/1/2等級)以解決此問題。
dirichlet也可在fit函數里使用,只要傳入pseudo_counts字典即可,如下面這樣

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) 
  • 1
  • 2

到此為止,模型已經完全構建完畢,下面可以開始使用其進行推理了。

推理

首先可以通過一些方法查看模型

#輸出節點信息 print(model.nodes()) #輸出依賴關系 print(model.edges()) #查看某節點概率分布 print(model.get_cpds('Pclass').values) 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

當然我們更關心的是給定某些節點后,感興趣節點的概率等,這就是推理。
貝葉斯網絡推理分成:
1.后驗概率問題:
表達為求P(Q|E=e) 其中Q為查詢變量 E為證據變量
即如本例子里已知一個人,女,<15歲,高票價,問生還幾率是多少?
2.最大后驗假設問題(MAP):
在這里插入圖片描述
已知證據E時,對某些變量的轉態組合感興趣(稱假設變量H),找最可能組合就是最大后驗假設問題。如本例子里,問一個活下來的女乘客最可能有是什么艙段,年齡?
3.最大可能解釋問題(MPE):是2的特例,即假設包含網絡里所有非證據變量(同時也可包含證據變量)

貝葉斯網絡推理主要有兩類方法,精確推理(變量化簡Variable Elination和置信傳播)和近似推理(如mcmc采樣),一般精確推理足以解決
pgmpy解決1可以用query函數 解決2,3可以用map_query函數
通過這些查詢可以獲得我們感興趣的關於因果關系信息,這是貝葉斯網絡模型的一大優勢。此處的因果關系並不可以解釋為一般意義上的邏輯因果,而是表示一種概率上的相關,比如我們不能將P(天亮了|公雞打鳴)很高解釋為是因為公雞打鳴天才亮的。

from pgmpy.inference import VariableElimination model_infer = VariableElimination(model) q = model_infer.query(variables=['Survived'], evidence={'Fare': 0}) print(q['Survived']) ''' +------------+-----------------+ | Survived | phi(Survived) | +============+=================+ | Survived_0 | 0.6341 | +------------+-----------------+ | Survived_1 | 0.3659 | +------------+-----------------+ ''' q = model_infer.map_query(variables=['Fare','Age','Sex','Pclass','Cabin'], evidence={'Survived': 1}) print(q)#{'Sex': 0, 'Fare': 0, 'Age': 1, 'Pclass': 2, 'Cabin': 0} 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16

上面的代碼使用了VariableElimination方法,亦可用BeliefPropagation,其有相同的接口。

與fit函數類似,也提供了輸入dataframe的簡便推理方法predict,如下
只要剔除想預測的列輸入predict函數,就可以得到預測結果的dataframe

predict_data=test.drop(columns=['Survived'],axis=1) y_pred = model.predict(predict_data) print((y_pred['Survived']==test['Survived']).sum()/len(test)) #測試集精度0.8131868131868132 
  • 1
  • 2
  • 3
  • 4

只是用了泰坦尼克數據集的一部分特征隨手設計網絡就可以達到不錯的效果了,精度81.3%,上傳kaggle的正確率0.77990。

自動設計網絡結構->使用結構學習方法

ref:
https://github.com/pgmpy/pgmpy_notebook
《貝葉斯網絡引論》

自動設計網絡結構的核心問題有兩個,一個是評價網絡好壞的指標,另一個是查找的方法。窮舉是不可取的,因為組合數太大,只能是利用各種啟發式方法或是限定搜索條件以減少搜索空間,因此產生兩大類方法,Score-based Structure Learning與constraint-based structure learning 以及他們的結合hybrid structure learning。
1.Score-based Structure Learning
這個方法依賴於評分函數,常用的有bdeu k2 bic,更合理的網絡評分更高,如下面的例子
此例子隨機產生x y 並令z=x+y,顯然X -> Z <- Y的結構合理

import pandas as pd import numpy as np from pgmpy.estimators import BdeuScore, K2Score, BicScore from pgmpy.models import BayesianModel # create random data sample with 3 variables, where Z is dependent on X, Y: data = pd.DataFrame(np.random.randint(0, 4, size=(5000, 2)), columns=list('XY')) data['Z'] = data['X'] + data['Y'] bdeu = BdeuScore(data, equivalent_sample_size=5) k2 = K2Score(data) bic = BicScore(data) model1 = BayesianModel([('X', 'Z'), ('Y', 'Z')]) # X -> Z <- Y model2 = BayesianModel([('X', 'Z'), ('X', 'Y')]) # Y <- X -> Z print(bdeu.score(model1)) print(k2.score(model1)) print(bic.score(model1)) print(bdeu.score(model2)) print(k2.score(model2)) print(bic.score(model2)) ''' -13936.101051153362 -14326.88012027081 -14292.1400887 -20902.744280734016 -20929.567083476162 -20946.7926535 ''' 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32

X -> Z <- Y 的評分更高
而依據評分函數進行搜索的搜索方法常用的有窮舉(5個節點以下可用)和 爬山算法(一個貪婪算法)pympy的實現如下:

from pgmpy.estimators import HillClimbSearch # create some data with dependencies data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH')) data['A'] += data['B'] + data['C'] data['H'] = data['G'] - data['A'] hc = HillClimbSearch(data, scoring_method=BicScore(data)) best_model = hc.estimate() print(best_model.edges()) #[('A', 'C'), ('A', 'B'), ('C', 'B'), ('G', 'A'), ('G', 'H'), ('H', 'A')] 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12

2.Constraint-based Structure Learning
比如根據獨立性得到最優結構的方法,相對來講前一種更有效

from pgmpy.independencies import Independencies ind = Independencies(['B', 'C'], ['A', ['B', 'C'], 'D']) ind = ind.closure() # required (!) for faithfulness model = ConstraintBasedEstimator.estimate_from_independencies("ABCD", ind) print(model.edges()) #[('A', 'D'), ('B', 'D'), ('C', 'D')] 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

回到泰坦尼克的例子,使用HillClimbSearch試了以下,在自己的測試集得到了和之前手工設計網絡相同的精度(kaggle測試成績略低一些),但是模型結構更復雜,不過看看給出的模型可以發現一些有趣的東西,比如Cabin Pclass Fare 相關,Age還影響了Sex等等

from pgmpy.estimators import HillClimbSearch from pgmpy.estimators import BdeuScore, K2Score, BicScore hc = HillClimbSearch(train, scoring_method=BicScore(train)) best_model = hc.estimate() #print(best_model.edges()) showBN(best_model) 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

在這里插入圖片描述

best_model.fit(train, estimator=BayesianEstimator, prior_type="BDeu") # default equivalent_sample_size=5 predict_data=test.drop(columns=['Survived'],axis=1) y_pred = best_model.predict(predict_data) (y_pred['Survived']==test['Survived']).sum()/len(test)#測試集精度 #0.8131868131868132 
  • 1
  • 2
  • 3
  • 4
  • 5

模型保存

先驗

這里摘錄一些下文對先驗的介紹
ref:Bayesian Methods for Hackers chp6

貝葉斯先驗可以分為兩類:客觀先驗和主觀先驗。客觀先驗的目的是讓數據對后驗影響最大,主觀先驗的目的是讓從業者對前驗表達自己的觀點。
事實上,從選擇先驗分布開始就已經開始搭建模型了,屬於建模的一部分。如果后驗不符合要求,自然可修改更換先驗,這是無關緊要的。沒有正確的模型,只有有用的模型。
經驗貝葉斯方法
一種叫經驗貝葉斯方法融合了頻率派的做法,采用 數據>先驗>數據>后驗 的方法,從觀測數據的統計特征構建先驗分布。這種方法實則違背了貝葉斯推理 先驗>數據>后驗 的思路。
從專家獲得先驗分布
可以考慮實驗轉盤賭法捕獲專家認為的先驗分布形狀。做法是讓專家吧固定總數的籌碼放在各個區間上以此來表達各個區段的概率,如下所示。之后可用合適的分布對此進行建模擬合,得到專家先驗。
在這里插入圖片描述
判斷先驗是否合適
只要先驗在某處概率不為零,后驗就有機會在此處表達任意的概率。當后驗分布的概率堆積在先驗分布的上下界時,那么很肯先驗是不大對的。(比如用Uniform(0,0.5)作為真實值p=0.7的先驗,那么后驗推斷的結果會堆積在0.5一側,表明真實值極可能大於0.5)

練手數據集

Binary Classification

Multiclass Classification

下載

鏈接: https://pan.baidu.com/s/1Im_Ls_Nul-swY9vrTueViw 提取碼: 7isb


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM