1. 例題1
藤田一弥. (2015). 見えないものをさぐる―それがベイズ ツールによる実践ベイズ統計. 株式会社 オーム社. P99 4-2 例題2を利用して、例題を見ていきます。以下のネットワークを考えます。W(Weather), R(Rain), S(Sprinkler), G(Ground)となっていて、天気の良し悪しで、地面が濡れているかどうかを、雨またはスプリンクラーの影響を入れてネットワークができています。
.png)
2. 例題2 : 涌井良幸, & 涌井貞美. (2012). 史上最強図解 これならわかる! ベイズ統計学. ナツメ社. P170 例題2を参照.png)
この図のようなネットワークを作成します。
藤田一弥. (2015). 見えないものをさぐる―それがベイズ ツールによる実践ベイズ統計. 株式会社 オーム社. P99 4-2 例題2を利用して、例題を見ていきます。以下のネットワークを考えます。W(Weather), R(Rain), S(Sprinkler), G(Ground)となっていて、天気の良し悪しで、地面が濡れているかどうかを、雨またはスプリンクラーの影響を入れてネットワークができています。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import networkx as nx | |
| import matplotlib.pyplot as plt | |
| # ベイジアンネットワークの有向グラフを作成 | |
| G = nx.DiGraph() | |
| # ノードを追加 | |
| G.add_node("W") | |
| G.add_node("S") | |
| G.add_node("R") | |
| G.add_node("G") | |
| # エッジを追加 | |
| G.add_edge("W", "S") | |
| G.add_edge("W", "R") | |
| G.add_edge("S", "G") | |
| G.add_edge("R", "G") | |
| # ベイジアンネットワークを描画 | |
| pos = nx.spring_layout(G) # グラフをレイアウト | |
| nx.draw(G, pos, with_labels=True, node_size=3000, node_color="skyblue", font_size=15, font_weight="bold") # グラフを描画 | |
| plt.show() |
.png)
この例題を通して実施したいことは、以下の3つです。
(1) ベイジアンネットワークのモデル(構造)を作る(ネットワークの形、尤度)
(2) 前方への確率(条件付き確率)を求める (P(G|W=w))など、グラフの左から右への確率
(3) 後方への確率(事後確率)※これがメインのもの、グラフの右から左への確率
今回は、pgmpyを使ってベイジアンネットワークでの分析をするので、インストールをしてください。
!pip install pgmpy
(1) ベイジアンネットワークのモデル(構造)を作る(ネットワークの形、尤度)
まず、例題で与えられたネットワークを作成していきます。また、推移確率も追加していきます。この推移確率を作成するときが間違えやすく、行がvariable(確率変数)、列がevidence(条件)の組み合わせになっているので、推移確率が今回の例題のように与えられた場合は、注意をしてください。今回の推移確率とWの事前分布はソースコードを確認してください。例えば、[[P(S=0|W=0), P(S=0|W=1)], [P(S=1|W=0), P(S=1|W=1)]]のように与えます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from pgmpy.models import BayesianNetwork | |
| from pgmpy.factors.discrete import TabularCPD | |
| # ベイジアンネットワークのモデルを作成 | |
| model = BayesianNetwork() | |
| # ノードを追加 | |
| model.add_nodes_from(['W', 'S', 'R', 'G']) | |
| # エッジを追加 | |
| model.add_edges_from([('W', 'S'), ('W', 'R'), ('S', 'G'), ('R', 'G')]) | |
| # 各ノードの条件付き確率分布を定義(行に変数、列に条件) | |
| cpd_w = TabularCPD(variable='W', variable_card=2, values=[[0.5], [0.5]]) #事前確率 | |
| cpd_s = TabularCPD(variable='S', variable_card=2, values=[[0.9, 0.5], [0.1, 0.5]], evidence=['W'], evidence_card=[2]) #[[P(S=0|W=0), P(S=0|W=1)], [P(S=1|W=0), P(S=1|W=1)]] | |
| cpd_r = TabularCPD(variable='R', variable_card=2, values=[[0.2, 0.8], [0.8, 0.2]], evidence=['W'], evidence_card=[2]) #[[P(R=0|W=0), P(R=0|W=1)], [P(R=1|W=0), P(R=1|W=1)]] | |
| cpd_g = TabularCPD(variable='G', variable_card=2, values=[[0.99, 0.1, 0.1, 0.01], [0.01, 0.9, 0.9, 0.99]], evidence=['S', 'R'], evidence_card=[2, 2]) #[[P(G=0|S=0, R=0), P(G=0|S=0, R=1), P(G=0|S=1, R=0), P(G=0|S=1, R=1)], [P(G=1|S=0, R=0), P(G=1|S=0, R=1), P(G=1|S=1, R=0), P(G=1|S=1, R=1)]] | |
| # モデルに条件付き確率分布を追加 | |
| model.add_cpds(cpd_w, cpd_s, cpd_r, cpd_g) | |
| # モデルの構造とCPDsが有効であるかをチェック | |
| assert model.check_model() | |
| # ベイジアンネットワークの概要を出力 | |
| print(model) |
(2) 前方への確率(条件付き確率)を求める (P(G|W=w))など、グラフの左から右への確率
今回は、「P(G = 1| W = 1)を計算：天気が晴れという条件のもとで、芝が濡れている確率」を求めていきます。つまり天気(W=1)が与えられたときのGの分布を求めていきます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #P101 (1) P(G = 1 | W = 1)を計算：天気が晴れという条件のもとで、芝が濡れている確率 | |
| from pgmpy.inference import VariableElimination | |
| # 推論エンジンの作成 | |
| infer = VariableElimination(model) | |
| # P(G = 1 | W = 1)を求める | |
| result = infer.query(variables=['G'], evidence={'W': 1}) | |
| print(result) | |
| probability_g_given_w_1 = result.values[1] # G = 1の確率を取得 | |
| print("P(G = 1 | W = 1):", probability_g_given_w_1) |
(3) 後方への確率(事後確率)※これがメインのもの、グラフの右から左への確率
今回は、「確率的推論(事後確率) P(W = 1 | G = 1) 芝が濡れているとき、晴れである確率」を求めていきます。「P(W = 1 | G = 1): 0.42538461538461536」が得られます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #P99 確率的推論(事後確率) P(W = 1 | G = 1) 芝が濡れているとき、晴れである確率は？ | |
| # P(W = 1 | G = 1)を求める | |
| result = infer.query(variables=['W'], evidence={'G': 1}) | |
| probability_w_given_g_1 = result.values[1] # W = 1の確率を取得 | |
| print("P(W = 1 | G = 1):", probability_w_given_g_1) |
もう一つ例題を見ていきます。
.png)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import networkx as nx | |
| import matplotlib.pyplot as plt | |
| # ベイジアンネットワークの有向グラフを作成 | |
| G2 = nx.DiGraph() | |
| # ノードを追加 | |
| G2.add_node("B") | |
| G2.add_node("E") | |
| G2.add_node("A") | |
| G2.add_node("P") | |
| G2.add_node("S") | |
| # エッジを追加 | |
| G2.add_edge("B", "A") | |
| G2.add_edge("E", "A") | |
| G2.add_edge("A", "P") | |
| G2.add_edge("A", "S") | |
| # ノードの位置を手動で設定 | |
| pos = {"B": (0.5, 2), "E": (1.5, 2), "A": (1, 1), "P": (0.5, 0), "S": (1.5, 0)} | |
| # ベイジアンネットワークを描画 | |
| #pos = nx.spring_layout(G2) # グラフをレイアウト | |
| nx.draw(G2, pos, with_labels=True, node_size=3000, node_color="skyblue", font_size=15, font_weight="bold") # グラフを描画 | |
| plt.show() |
このとき、B(Burglar)、E(EarthQuake)、A(Alert)、P(Police)、S(Security)といったノードになっています。泥棒(B)が入ったとき、アラート(A)がなり、セキュリティ会社(S)に連絡がいくという流れです。
(1) ベイジアンネットワークのモデル(構造)を作る(ネットワークの形、尤度)
これは例題で与えられているので、尤度(推移確率)の構造に注意して作成していきます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from pgmpy.models import BayesianNetwork | |
| from pgmpy.factors.discrete import TabularCPD | |
| # ベイジアンネットワークのモデルを作成 | |
| model = BayesianNetwork() | |
| # ノードを追加 | |
| model.add_nodes_from(['B', 'E', 'A', 'P', 'S']) | |
| # エッジを追加 | |
| model.add_edges_from([('B', 'A'), ('E', 'A'), ('A', 'P'), ('A', 'S')]) | |
| # 各ノードの条件付き確率分布を定義 | |
| cpd_b = TabularCPD(variable='B', variable_card=2, values=[[0.99], [0.01]]) # Bの事前確率 | |
| cpd_e = TabularCPD(variable='E', variable_card=2, values=[[0.98], [0.02]]) # Eの事前確率 | |
| cpd_a = TabularCPD(variable='A', variable_card=2, | |
| values=[[0.92, 0.74, 0.06, 0.05], # P(A|B, E) | |
| [0.08, 0.26, 0.94, 0.95]], | |
| evidence=['B', 'E'], evidence_card=[2, 2]) | |
| cpd_p = TabularCPD(variable='P', variable_card=2, values=[[0.5, 0.5], [0.5, 0.5]], evidence=['A'], evidence_card=[2]) # Pの条件付き確率(今回は未設定) | |
| cpd_s = TabularCPD(variable='S', variable_card=2, values=[[0.9, 0.3], [0.1, 0.7]], evidence=['A'], evidence_card=[2]) # Sの条件付き確率 | |
| # モデルに条件付き確率分布を追加 | |
| model.add_cpds(cpd_b, cpd_e, cpd_a, cpd_p, cpd_s) | |
| # モデルの構造とCPDsが有効であるかをチェック | |
| assert model.check_model() | |
| # ベイジアンネットワークの概要を出力 | |
| print(model) |
(2)前方への確率(条件付き確率)を求める (P(G|W=w))など、グラフの左から右への確率
これは演習としますので、やってみてください。
(3) 後方への確率(事後確率)※これがメインのもの、グラフの右から左への確率
「P(B | S) セキュリティ会社に連絡がいったとき、泥棒が入った確率」を求めてみます。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #P(B | S) セキュリティ会社に連絡がいったとき、泥棒が入った確率 | |
| from pgmpy.inference import VariableElimination | |
| # 推論エンジンの作成 | |
| infer = VariableElimination(model) | |
| # P(B | S)を求める | |
| result = infer.query(variables=['B'], evidence={'S': 1}) | |
| probability_b_given_s = result.values[1] # B = 1の確率を取得 | |
| print("P(B = 1 | S = 1):", probability_b_given_s) |
