ベイジアンネットワーク

ベイジアンネットワーク#

ベイジアンネットワーク（Bayesian network, Bayesnet）は複数の確率変数の間の依存関係をDAGで表現するモデル

分析の流れ#

確率変数間の関係性の表現#

確率変数間の関係を確率の積で表す。例えばX→Yという依存関係（因果関係）は

P(X,Y)=P(Y|X)P(X)

と分解されると考える。P(Y|X)はXの値によって決まるので、Xの値に応じてYが決まるという依存関係を表す。

../_images/1768e8890c91d214cb4531762e1e9848af3e8a609b92a209d715490b98b31f06.svg

上の図の場合、

P(Y,X1,X2)=P(Y|X1,X2)P(X1)P(X2)

となる

この条件付き確率を確率変数間の関係性として扱い、 条件付き確率表 （Conditional Probabilities Tables: CPT）を使って表現する。

（ CPD (Conditional Probability Distribution) という呼び方もある様子？）

CPTは単に確率変数の値と確率の表である。例えば以下のようになる

x1 P(x1)
0 0 0.6
1 1 0.4
x2 P(x2)
0 0 0.4
1 1 0.6
Hide code cell source 
import numpy as np
import pandas as pd

cpt_y = pd.DataFrame([
    {"x1": 0, "x2": 0, "y": 0, "P(y=0|x1,x2)": 0.2},
    {"x1": 0, "x2": 1, "y": 0, "P(y=0|x1,x2)": 0.3},
    {"x1": 1, "x2": 0, "y": 0, "P(y=0|x1,x2)": 0.4},
    {"x1": 1, "x2": 1, "y": 0, "P(y=0|x1,x2)": 0.1},
])
del cpt_y["y"]
cpt_y["P(y=1|x1,x2)"] = 1 - cpt_y["P(y=0|x1,x2)"]
cpt_y
x1 x2 P(y=0|x1,x2) P(y=1|x1,x2)
0 0 0 0.2 0.8
1 0 1 0.3 0.7
2 1 0 0.4 0.6
3 1 1 0.1 0.9

パラメータの推定#

最尤推定でCPTの各セルにいれるべき確率（これをパラメータθと表すことにする）を推定する。

確率変数のインデックスをi、確率変数の値をj、条件付き確率については条件をkとすると、推定したいパラメータはθi,j,kとなる

例えばP(x=1)はx=1を満たすサンプル数N1と全体のサンプル数の比率

P(x=1)=N1N

で求められるように、

同様に、条件i,j,kを満たすサンプル数をNi,j,kと書くことにすると

θi,j,k=Ni,j,kN

となる

モデルの評価#

データに対するネットワークの当てはまりの良さの指標にはAICやBICなどが使われる

BICは

BIC=−2l(θ|X)+k(log⁡N)

Xはデータ、kはモデルのパラメータ数、Nはサンプルサイズ、l(θ|X)はモデルの対数尤度

分析例（スクラッチ）#

データの用意#

まずデータを生成する。上記のCPTを真の値として、そこから疑似乱数でデータを生成する。

import numpy as np

n = 100
np.random.seed(0)

x1 = np.random.binomial(n=1, p=cpt_x1.query("x1 == 1")["P(x1)"].iloc[0], size=n)
x2 = np.random.binomial(n=1, p=cpt_x2.query("x2 == 1")["P(x2)"].iloc[0], size=n)
df = pd.DataFrame({"x1": x1, "x2": x2})

# yを生成
for x1 in [0, 1]:
    for x2 in [0, 1]:
        q = f"x1 == {x1} & x2 == {x2}"
        values = np.random.binomial(
            n=1,
            p = cpt_y.query(q)["P(y=1|x1,x2)"].iloc[0],
            size = df.query(q).shape[0],
        )
        df.loc[df.query(q).index, "y"] = values
df = df.astype({"y": "int"})

df.tail(3)
x1 x2 y
97 0 1 1
98 1 1 0
99 0 1 1

パラメータの推定#

パラメータはθi,j,kを求めていく

θi,j,k=Ni,j,kN

確率変数のインデックスをi、確率変数の値をj、条件付き確率の条件がkである。

P(x1=0)については、条件付き確率ではないのでk=[]と表すことにするとθ1,0,[]となる

J = [0, 1]

# θ_{1,j,[]}
theta_x1 = []
# サンプル数 N_ijkも保存する
n_x1 = []

for j in J:
    n_ij = df.query(f"x1 == {j}").shape[0]
    theta_x1.append(n_ij / n)
    n_x1.append(n_ij)

theta_x1

[0.62, 0.38]

# θ_{2,j,[]}
theta_x2 = []
# サンプル数 N_ijkも保存する
n_x2 = []

for j in J:
    n_ij = df.query(f"x2 == {j}").shape[0]
    theta_x2.append(n_ij / n)
    n_x2.append(n_ij)

theta_x2

[0.44, 0.56]

# 条件をリストアップ
k_values = []
for x1 in [0, 1]:
    for x2 in [0, 1]:
        k_values.append((x1, x2))
k_values

[(0, 0), (0, 1), (1, 0), (1, 1)]

# θ_{3,j,[]}
theta_y = np.zeros(shape=(2, 4))  # j * k
n_y = np.zeros(shape=(2, 4))  # j * k
for k, (x1, x2) in enumerate(k_values):
    q = f"x1 == {x1} & x2 == {x2}"
    for j in J:
        n_ijk = df.query(q).query(f"y == {j}").shape[0]
        theta_y[j, k] = n_ijk / n
        n_y[j, k] = n_ijk
theta_y

array([[0.06, 0.08, 0.08, 0.03],
       [0.23, 0.25, 0.07, 0.2 ]])

モデルの評価#

モデルの対数尤度l(θ|X)に比例する値を計算する

l(θ∣X)∝∑i∑j∑k(Nijk)log⁡θ^i,j,k
thetas = [theta_x1, theta_x2, theta_y]
nums = [n_x1, n_x2, n_y]

log_likelihood = 0
for i in range(len(thetas)):
    for j in [0, 1]:
        if i == 2:
            for k in range(len(k_values)):
                log_likelihood += nums[i][j][k] * np.log(thetas[i][j][k])
        else:
            log_likelihood += nums[i][j] * np.log(thetas[i][j])
log_likelihood

-322.07467382966473
BIC=−2l(θ|X)+k(log⁡N)

Xはデータ、kはモデルのパラメータ数、Nはサンプルサイズ、l(θ|X)はモデルの対数尤度

# 『Pythonによる因果分析』で見てもkの値がよくわからなかった
k = len(thetas) * len(J)
bic = -2 * log_likelihood + k * np.log(n)
bic

671.780368775258

pympyパッケージによる推定#

Bayesian Network — pgmpy 0.1.23 documentation

from pgmpy.models import BayesianNetwork

DAG = [('x1', 'y'), ('x2', 'y')] # x1 -> y, x2 -> y
model = BayesianNetwork(DAG)

model.fit(df)
model.get_cpds()

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[14], line 4
      1 from pgmpy.models import BayesianNetwork
      3 DAG = [('x1', 'y'), ('x2', 'y')] # x1 -> y, x2 -> y
----> 4 model = BayesianNetwork(DAG)
      6 model.fit(df)
      7 model.get_cpds()

File /usr/local/lib/python3.10/site-packages/pgmpy/models/BayesianNetwork.py:3, in BayesianNetwork.__init__(self, ebunch, latents, lavaan_str, dagitty_str)
      2 def __init__(self, ebunch=None, latents=set(), lavaan_str=None, dagitty_str=None):
----> 3     raise ImportError(
      4         "BayesianNetwork has been deprecated. Please use DiscreteBayesianNetwork instead."
      5     )

ImportError: BayesianNetwork has been deprecated. Please use DiscreteBayesianNetwork instead.

# 推定されたCPT
print(model.get_cpds()[0])

+-------+------+
| x1(0) | 0.62 |
+-------+------+
| x1(1) | 0.38 |
+-------+------+

from pgmpy.estimators import BicScore
bic = BicScore(df)
print(bic. score(model))

-201.14661436972517

参考文献#

Contents