PCアルゴリズム

PCアルゴリズム#

Peter-Clark (PC) algorithm

Step1:すべてのノード間が無向のエッジで結ばれたグラフを作る

../../_images/66d74c62e0ee23f874ee9f80a4ca59210ca98972766ceedc6b508fbca3dd8507.svg

Step 2:条件付き独立性を検定し、不要なエッジの削除

独立性や条件付き独立性をカイ二乗検定で判定していく。 \(X_i \perp X_j\) または \(X_i \perp X_j \mid S\) が統計的に判定できた場合、そのエッジを削除し、分離集合(separating set)\(S\) を記録する。

例:\(X_1 \perp X_3 \mid X_2\) が成立 → 辺 X1–X3 を削除、sep(X1,X3) = {X2}

../../_images/43d653a47ff27f80c027a1ad7f9474fcc968f14a9a393b5791f76b06845ccc79.svg

Step 3:V-構造(Collider)の同定

3変数 \(X_i - X_j - X_k\) のパターンで、

  • \(X_i\)\(X_k\) が「隣接していない」(エッジがない、\(X_i \perp X_k \mid S\)

  • その分離集合 sep(\(X_i, X_k\)) に \(X_j\) が含まれない(\(X_j \notin S\)

のとき、\(X_i \rightarrow X_j \leftarrow X_k\)(V構造、collider) と向きを付ける。

../../_images/ca5d4bc321796f96bea0e5b1992b9b35abbccdd779943213853fe1e3b9549656.svg

Step 4:向き付け規則(Orientation Rules)で残りの辺を決める

追加のルールを適用し、サイクルを作らない・新たな collider を作らないように残りの無向辺に向きを付ける。

代表的ルール:

  • Orientation propagation:\(X \to Y - Z\) なら \(Y \to Z\)

  • Acyclicity constraint(サイクル禁止)

  • Avoid new colliders(既存以外の新規 collider を作らない)

../../_images/9980a157f07014848566e1a98650c4bd93e190c0c75f1b540e654894e204ffca.svg

限界・注意点#

PCアルゴリズムは以下の点に注意が必要

1️⃣すべての変数間の因果の方向が定まらないことがある

同じ条件付き独立性を与える因果グラフの集合を マルコフ同値類(Markov equivalence class: MEC) という。
PCアルゴリズムではMECの中から因果グラフを一意に特定できない場合がある。

2️⃣サンプルサイズが必要

サンプルサイズが小さい場合はカイ二乗分布を用いた近似精度が悪化するため、適切に検定できなくなる

3️⃣高次元データに弱い

PCアルゴリズムは条件付き独立性の検定を実施する順番に影響されやすい。そのため次元数が多いと推定が不安定になりやすい。 これを改善したPC-stableというアルゴリズムもある。

4️⃣未観測の交絡因子と選択バイアスが考慮されていない

因果的十分性の仮定をおいているため、この仮定が満たされない場合は因果探索の推定精度が下がる。
この点を改良したFCIアルゴリズムが提案されている。

5️⃣検定の有意水準の設定によって推定結果が変わる

6️⃣計算量が多い

すべての変数の組み合わせに対して統計的検定を行うため、変数の数が多いと処理しきれない

実装#

pgmpyパッケージに実装されている

PC (Constraint-Based Estimator) — 1.0.0 | pgmpy docs

# データセットの生成
import numpy as np
from pgmpy.utils import get_example_model
np.random.seed(42)
model = get_example_model("asia")
model.seed = 42
df = model.simulate(int(1e3))

# 因果構造を推定
from pgmpy.estimators import PC
dag = PC(data=df).estimate(ci_test='chi_square', return_type='dag')

# daftパッケージで作図
dag.to_daft().render()
INFO:numexpr.utils:NumExpr defaulting to 4 threads.
INFO:pgmpy: Datatype (N=numerical, C=Categorical Unordered, O=Categorical Ordered) inferred from data: 
 {'lung': 'C', 'tub': 'C', 'asia': 'C', 'xray': 'C', 'either': 'C', 'smoke': 'C', 'bronc': 'C', 'dysp': 'C'}
<Axes: >
../../_images/f5dcce4a59ce21334d97deecdb0daa42da2c684dd17fc7fe5473c593cdbbc0fe.png

参考#