Mímisbrunnr知恵の泉

← 因果推論 一覧

🎓 レベル：発展　|　重要度：C（発展的）　|　⚠️ この領域は急速に動くため随所に「要最新確認」と明示する

📎 前提：因果ダイアグラムとd分離　|　バックドア基準と識別　|　数理：構造方程式モデル・パス解析（統計）

要点（BLUF）

1. 因果探索とは——矢印を引く逆問題

第1〜7章では DAG を分野知識で与え、そこからバックドア基準で調整集合を選んだ（バックドア基準と識別）。因果探索はその逆で、観測データだけから矢印の向きを推定しようとする。「相関のパターン（条件付き独立の集合）」を手がかりに、それを生み得る DAG を探す。

これは本質的に難しい逆問題である。なぜなら因果ダイアグラムとd分離で見たとおり、異なる DAG が同じ条件付き独立を生むことがあるからだ。鎖 X→Z→YX\to Z\to Y と分岐 X←Z→YX\leftarrow Z\to Y はどちらも「X⊥Y∣ZX\perp Y\mid Z、X⊥̸YX\not\perp Y」を生み、データからは区別できない。

flowchart LR
  subgraph eq["マルコフ同値（データから区別不能）"]
    direction LR
    A1["X"] --> B1["Z"] --> C1["Y"]
    B2["Z"] --> A2["X"]
    B2 --> C2["Y"]
  end
  subgraph vs["v字構造（向きが識別できる）"]
    A3["X"] --> B3["Z"]
    C3["Y"] --> B3
  end

左の鎖と分岐は マルコフ同値——同じ独立関係を持つので観測データだけでは向きが決まらない。右の v字構造（コライダー） X→Z←YX\to Z\leftarrow Y だけは「X⊥YX\perp Y なのに X⊥̸Y∣ZX\not\perp Y\mid Z」という他にない指紋を持ち、向きが識別できる。因果探索の限界と可能性はすべてこの違いに集約される。

2. 4系統の概観と、それぞれが要る仮定

系統代表アルゴリズム中心となる仮定出力
制約ベースPC・FCI因果マルコフ性＋忠実性マルコフ同値類（CPDAG）
スコアベースGES因果マルコフ性＋忠実性＋スコア（BIC等）マルコフ同値類（CPDAG）
関数ベースLiNGAM線形＋非ガウスノイズ＋非巡回＋未観測交絡なし完全な DAG（全向き）
連続最適化NOTEARS線形 SEM＋スパース性＋微分可能な非巡回制約重み付き隣接行列（DAG）

3. 識別の仮定——マルコフ性・忠実性・非ガウス性

因果探索が依拠する基本仮定は 2 つ。

  1. 因果マルコフ性：各変数は親で条件づければ非子孫と独立。これにより同時分布が
P(V1,…,Vd)=∏i=1dP ⁣(Vi∣pa(Vi))P(V_1,\dots,V_d) = \prod_{i=1}^{d} P\!\left(V_i \mid \mathrm{pa}(V_i)\right)

と分解される（因果ダイアグラムとd分離）。

  1. 忠実性（faithfulness）：データに現れる条件付き独立は、すべて DAG の d 分離から来る（偶然の係数相殺で生じた独立はない）。これにより「独立 ⇔ d分離」が双方向に結ばれ、独立検定の結果から構造を逆算できる。

この 2 つだけからは、前節のとおり マルコフ同値類 までしか分からない。ここが識別の壁である。壁を越えるには分布の形に追加の仮定を入れる。LiNGAM の非ガウス性がその代表で、線形非ガウス非巡回モデル

x=Bx+e,x=(I−B)−1ex = Bx + e, \qquad x = (I-B)^{-1} e

（BB は下三角に並べ替え可能な係数行列、ee は互いに独立で非ガウスな外生ノイズ）を仮定すると、DAG が一意に識別できる（Shimizu et al. 2006）。ガウスだと回転不変性で向きが消えるが、非ガウスだと独立成分分析（ICA）の一意性が効いて向きが復元する。

⚠️ 要最新確認：因果探索はこの数年で連続最適化（NOTEARS 系）・スコアと制約の融合・未観測交絡対応（FCI/LiNA）・時系列版（VAR-LiNGAM）と急速に展開している。下記の手実装は「核となる識別の論理」を見せるためのもので、実務では後述ライブラリの最新版を使い、仮定の成否を必ず吟味すること。

4. コード①：PC 手続きの核——骨格・コライダー・同値類の壁

PC の中核だけを 3 変数で手実装する。条件付き独立検定は偏相関の Fisher z 検定で行う。同じ手続きを (A) コライダー X→Z←YX\to Z\leftarrow Y から生成したデータ、(B) 鎖 X→Z→YX\to Z\to Y から生成したデータに当て、コライダーは向きが決まる／鎖は同値類のままを確かめる。

# === PC手続き：条件付き独立検定で骨格→コライダー向き付け ===
import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
n = 4000


def partial_corr_test(data, i, j, cond):
    # data[:,i] と data[:,j] を cond（列番号リスト）で条件づけた偏相関のFisher z検定
    x = data[:, i]
    y = data[:, j]
    if len(cond) == 0:
        r = np.corrcoef(x, y)[0, 1]
    else:
        G = np.column_stack([np.ones(len(x))] + [data[:, k] for k in cond])
        bx, _, _, _ = np.linalg.lstsq(G, x, rcond=None)
        by, _, _, _ = np.linalg.lstsq(G, y, rcond=None)
        rx = x - G @ bx
        ry = y - G @ by
        r = np.corrcoef(rx, ry)[0, 1]
    dof = len(x) - len(cond) - 3
    r = np.clip(r, -0.9999, 0.9999)
    z = 0.5 * np.log((1 + r) / (1 - r)) * np.sqrt(dof)
    p = 2 * (1 - stats.norm.cdf(abs(z)))
    return r, p


def pc_three_nodes(data, names, alpha=0.01):
    # 3ノード専用の簡易PC：骨格推定→sepset記録→コライダー向き付け
    nodes = [0, 1, 2]
    adj = {(0, 1): True, (0, 2): True, (1, 2): True}
    sepset = {}
    for (i, j) in list(adj.keys()):
        others = [k for k in nodes if k != i and k != j]
        for cond in [[], others]:            # 空集合 と もう1つのノードで条件づけ
            r, p = partial_corr_test(data, i, j, cond)
            if p > alpha:                    # 独立と判断 → 辺を消し sepset を記録
                adj[(i, j)] = False
                sepset[(i, j)] = cond
                break
    arrows = []                              # コライダー判定：sepsetに無い共通隣接へ i→k←j
    for (i, j) in [(0, 1), (0, 2), (1, 2)]:
        if not adj[(i, j)]:
            k = [x for x in nodes if x != i and x != j][0]
            if adj.get((min(i, k), max(i, k)), False) and adj.get((min(j, k), max(j, k)), False):
                if k not in sepset.get((i, j), []):
                    arrows.append((i, k))
                    arrows.append((j, k))
    return adj, arrows


def describe(adj, arrows, names):
    edges = [f"{names[i]}-{names[j]}" for (i, j) in [(0, 1), (0, 2), (1, 2)] if adj[(i, j)]]
    dirs = [f"{names[a]}->{names[b]}" for (a, b) in arrows]
    return "骨格(無向辺): " + ", ".join(edges) + " / 向き付け: " + (", ".join(dirs) if dirs else "なし(同値類のまま)")


names = ["X", "Y", "Z"]

# --- (A) コライダー X->Z<-Y からデータ生成 ---
X = rng.normal(0, 1, n)
Y = rng.normal(0, 1, n)
Z = 1.0 * X + 1.0 * Y + rng.normal(0, 1, n)
adj, arrows = pc_three_nodes(np.column_stack([X, Y, Z]), names)
print("[コライダー X->Z<-Y で生成]")
print("  ", describe(adj, arrows, names))

# --- (B) 鎖 X->Z->Y からデータ生成 ---
X2 = rng.normal(0, 1, n)
Z2 = 1.0 * X2 + rng.normal(0, 1, n)
Y2 = 1.0 * Z2 + rng.normal(0, 1, n)
adj2, arrows2 = pc_three_nodes(np.column_stack([X2, Y2, Z2]), names)
print("[鎖 X->Z->Y で生成]")
print("  ", describe(adj2, arrows2, names))

出力：

[コライダー X->Z<-Y で生成]
   骨格(無向辺): X-Z, Y-Z / 向き付け: X->Z, Y->Z
[鎖 X->Z->Y で生成]
   骨格(無向辺): X-Z, Y-Z / 向き付け: なし(同値類のまま)

出力の意味：両者とも骨格は同じ「X–Z, Y–Z」（X,YX,Y は非隣接）。違いは向き付けに出る。

5. コード②：非ガウス性で同値類を脱出する（LiNGAM の核）

前節の壁を越える鍵が非ガウス性。LiNGAM の発想は単純で、真の向きでは「原因」と「回帰の残差」が独立になるが、逆向きでは独立にならない（非ガウスのとき）。独立度をカーネル統計量 HSIC（小さいほど独立）で測り、HSIC(原因, 残差) が小さい方を原因→結果と判定する。同じ手続きをガウス・ノイズに当てると両向きがほぼ同点になり、識別できないことも確かめる。

# === 非ガウス性で向きを決める（pairwise LiNGAM の核） ===
import numpy as np

rng = np.random.default_rng(3)
n = 1500


def hsic(a, b):
    # ガウスカーネルHSIC：a,b の独立性が高いほど0に近い（中央値ヒューリスティックで帯域決定）
    a = (a - a.mean()) / a.std()
    b = (b - b.mean()) / b.std()
    Da = (a[:, None] - a[None, :]) ** 2
    Db = (b[:, None] - b[None, :]) ** 2
    sa = np.median(Da[Da > 0])
    sb = np.median(Db[Db > 0])
    K = np.exp(-Da / sa)
    L = np.exp(-Db / sb)
    m = len(a)
    H = np.eye(m) - np.ones((m, m)) / m
    return float(np.trace(K @ H @ L @ H)) / (m - 1) ** 2


def direction_score(cause, effect):
    # cause で effect を線形回帰した残差と cause の従属度（HSIC）。真の向きほど小さい
    beta = np.cov(cause, effect)[0, 1] / np.var(cause)
    resid = effect - beta * cause
    return hsic(cause, resid)


def decide(u, v):
    s_uv = direction_score(u, v)   # u->v 仮説
    s_vu = direction_score(v, u)   # v->u 仮説
    winner = "X->Y" if s_uv < s_vu else "Y->X"
    return s_uv, s_vu, winner


# --- (A) 非ガウス（一様）ノイズ：真の構造は X->Y ---
X = rng.uniform(-1, 1, n)
Y = 0.8 * X + rng.uniform(-1, 1, n)
s_uv, s_vu, win = decide(X, Y)
print("[非ガウス・ノイズ / 真の向き X->Y]")
print("  HSIC(X->Y仮説):", round(s_uv, 5), " HSIC(Y->X仮説):", round(s_vu, 5), " 判定:", win)

# --- (B) ガウス・ノイズ：真の構造は X->Y だが向きは識別不能のはず ---
Xg = rng.normal(0, 1, n)
Yg = 0.8 * Xg + rng.normal(0, 1, n)
s_uv2, s_vu2, win2 = decide(Xg, Yg)
ratio = max(s_uv2, s_vu2) / min(s_uv2, s_vu2)
print("[ガウス・ノイズ / 真の向き X->Y だが識別困難]")
print("  HSIC(X->Y仮説):", round(s_uv2, 5), " HSIC(Y->X仮説):", round(s_vu2, 5), " 判定:", win2)
print("  ガウスでの両向きHSIC比(1に近いほど見分け不能):", round(ratio, 3))

出力：

[非ガウス・ノイズ / 真の向き X->Y]
  HSIC(X->Y仮説): 8e-05  HSIC(Y->X仮説): 0.00798  判定: X->Y
[ガウス・ノイズ / 真の向き X->Y だが識別困難]
  HSIC(X->Y仮説): 0.00012  HSIC(Y->X仮説): 0.0001  判定: Y->X
  ガウスでの両向きHSIC比(1に近いほど見分け不能): 1.191

出力の意味：

6. ライブラリ（要最新確認）

実務では手実装でなく専用ライブラリを使う。いずれも API が動きやすいので必ず最新ドキュメントを確認すること。本ノートは上の手実装だけで完結しており、ライブラリは補助。

⚠️ 導入時の注意：これらは古い numpy/scipy/pandas に固定しがちで、最新環境ではダウングレードや依存衝突を起こすことがある。隔離環境（仮想環境）に入れ、まず本ノートの手実装で核を理解してから使うと安全。

7. 直観——なぜ非ガウスだと向きが見えるのか

ガウス分布は 2 次（共分散）までで決まり、X→YX\to Y と Y→XY\to X は同じ共分散構造を生むので回帰の残差はどちらでも原因と無相関かつ独立になる。ところが非ガウスだと、残差が原因と独立になるのは真の向きだけ。逆向きの残差は原因と「無相関だが従属」になり、その従属が高次モーメント（非ガウス性）に現れる。HSIC はその従属を検出している。つまり非ガウス性は向きの情報を分布の形に刻む。LiNGAM はこの刻印を読み、ICA の一意性で全 DAG を復元する。

⚠️ よくある誤解・落とし穴

関連ノート