ベイジアンネットワーク入門:pgmpyによる因果探索と因果推論の実践

近年の機械学習(ML)モデルは、ビッグデータ解析において非常に高い予測精度を達成しています。しかし、その意思決定に至るプロセスが不透明な「ブラックボックス」となってしまう課題が指摘されています。データから相関関係を発見するだけでは、「なぜその結果が導かれたのか」という因果関係の明確な理解には至りません。

今後、AIシステムに求められる透明性(Transparency)、公平性(Fairness)、そして頑健性(Robustness)といった重要な要素を確保するためには、観測データから因果効果を識別し、その構造を理解することが不可欠です。

今回解説するベイジアンネットワーク(Bayesian Network, BN)は、この因果関係をモデル化するための理論的基盤を提供します。BNは、変数間の依存関係を有向非巡回グラフ(DAG)と条件付き確率で効率的に表現する確率モデルです。

そして、この理論をPython上で柔軟かつ包括的に実践するためのライブラリが、pgmpy(Probabilistic Graphical Models in Python)です。pgmpyは、モデル構造の学習(因果探索)や、do演算子を用いた因果効果の識別・推定(因果推論)など、幅広い機能をサポートしています。

ベイジアンネットワークと因果推論の基礎

ベイジアンネットワークは、確率的グラフィカルモデル(Probabilistic Graphical Model, PGM)の一種であり、複雑なシステムを構成する変数間の確率的な依存関係をグラフィカルに、かつ効率的に表現するためのフレームワークです。

このベイジアンネットワークが、単なる確率モデルに留まらず、因果推論(Causal Inference)の文脈で重要視されるのは、そのグラフ構造が因果構造をモデル化するための厳密な基盤を提供するためです。

ベイジアンネットワークの構造と数理的基礎

ベイジアンネットワークの構造は、有向非巡回グラフ(Directed Acyclic Graph, DAG)によって定義されます。

  • ノード (Nodes):ネットワークを構成する各ノードは、分析対象となるシステムの確率変数を表します。例えば、「天気」や「病気の有無」などがノードに対応します。
  • エッジ (Edges):ノード間を結ぶ有向矢印(エッジ)は、変数間の直接的な確率的依存関係を示します。因果的ベイジアンネットワークにおいては、この矢印は「原因から結果へ」の一方向的な因果関係を表すものと解釈されます。
  • 非巡回性 (Acyclicity):DAGには、矢印に沿って進んだときに元のノードに戻ってくる閉路(サイクル)が存在しないという制約があります。これは、ある事象が自分自身の原因とはならないという因果関係の基本的な前提を反映しています。

このグラフ構造により、システム内の全変数の同時確率分布を、各ノードの条件付き確率(Conditional Probability, CPT)の積として効率的に分解表現できることが、BNの最大の強みの一つです。

$$\mathbb{P}(X_{1},…,X_{n}) = \prod_{i=1}^{n}\mathbb{P}(X_{i} | \text{Parents}(X_{i})) \tag{1}$$

この分解は、同時確率分布のパラメータ数を大幅に削減し、計算論的・統計的に扱いやすいモデルを実現します。

グラフ構造から独立性を読み解く:d-分離

BNのグラフ構造が確立されると、変数間の条件付き独立性(Conditional Independence)の関係を、統計的な計算を介さずに、グラフィカルな基準で読み解くことができます。これがd-分離(有向分離)です。

d-分離は、特定のノード集合 \(Z\) を条件としたときに、他の2つのノード集合 \(X\) と \(Y\) の間のすべてのパスが「ブロック」されているかどうかを判定するルールです。パスがブロックされている場合、 \(X\) と \(Y\) は \(Z\) のもとで条件付き独立である (\(X \perp Y | Z\)) と結論づけられます。

d-分離を理解する上で重要なのは、パスを構成する基本的な3つの構造の振る舞いです。

構造パスの振る舞い因果推論における役割
チェーン (Chain)\(A \rightarrow B \rightarrow C\)中間変数 \(B\) を条件付けない限り「オープン」。 \(B\) を調整すると \(A\) と \(C\) は独立になる(媒介変数)。
フォーク (Fork)\(A \leftarrow B \rightarrow C\)共通原因 \(B\)を条件付けない限り「オープン」。\(B\) を調整すると \(A\) と \(C\) は独立になる(交絡因子の調整)。
コライダー (Collider)\(A \rightarrow B \leftarrow C\)共通結果 \(B\) を条件付けない限り「ブロック」されている。\(B\) を条件付けるとパスが「オープン」になり、見せかけの相関が生じる(選択バイアス)。

このd-分離の概念は、因果推論における交絡(Confounding)や選択バイアス(Selection Bias)といった核心的な問題を、構造的に特定し、分析において不適切な変数を調整してしまう過ちを避ける上で極めて重要です。

因果探索と因果推論

因果分析は大きく分けて、因果構造をデータから見つける因果探索と、その構造を用いて介入効果を推定する因果推論に分けられます。

1. 因果探索 (Causal Discovery)

因果探索は、大量の観測データから、変数間に存在する因果構造 (DAG) そのものを学習・推定するプロセスです。ドメイン知識が不完全な場合に特に重要となります。

この構造学習には主に二つの手法が用いられます。

  • 制約ベース手法 (Constraint-Based):データが持つ統計的な条件付き独立性のテスト結果を利用して、DAGの骨格やエッジの向きを決定します(例:PCアルゴリズム)。
  • スコアベース手法 (Score-Based):候補となるDAG構造に対してスコア関数(例:BIC, BDeu)を計算し、そのスコアを最大化するようにグラフ空間を探索します(例:GES, Hill-Climb Search)。最近では、ベイズ最適化を用いて有望なDAG構造を効率的に発見する手法(例:DrBO)も提案され、計算効率と精度の向上が図られています。

2. 因果推論 (Causal Inference)

因果推論が目指すのは、「もし我々が \(X\) の値を強制的に \(X\) に設定したら、\(Y\) はどうなるか」という介入(Intervention)の効果を定量化することです。これは、単に「\(X\) が \(x\) であると観測されたときの \(Y\) の確率」\(P(Y|X=x)\) を求める従来の予測とは根本的に異なります

  • \(do\) 演算子の導入:Judea Pearlは、この「介入」を数学的に表現するために\(do\) 演算子を導入し、因果効果を \(P(Y|\text{do}(X=x))\) と表記しました。
  • グラフィカルな表現:\(do(X=x)\) という介入は、因果グラフ上で「グラフ手術」として表現され、\(X\) に向かう全ての矢印を削除し、\(X\) の値を \(x\) に固定する操作に対応します。これにより、\(X\) がその通常の原因から切り離され、外部からの強制的な介入であることがモデル化されます。
  • 因果効果の識別基準:観測データから因果効果を識別するための最も実用的で広く使われる基準がバックドア基準(Backdoor Criterion)です。これは、特定の変数集合 $Z$ が $X$ から $Y$ への因果効果を識別するための条件を満たす場合に、バックドア調整公式を用いて交絡の影響を取り除き、純粋な因果効果を推定するものです。

$$P(Y=y|\text{do}(X=x))=\sum_{z}P(Y=y|X=x,Z=z)P(Z=z) \tag{2}$$

実用例:因果モデルの価値

因果モデルの構築は、予測精度向上だけでなく、AIシステムの透明性信頼性を高めます。特に、複雑な要因が絡み合う実世界の問題において、その価値を発揮します。

  • 医療・ヘルスケア:臨床データやゲノムデータなどの異種混合データを統合し、疾患の診断支援や予後予測、治療効果の推定モデル開発に利用されます。グラフ構造は推論の根拠を視覚的に示し、「ホワイトボックス」なAIとして機能します。
  • 社会科学・公平性:アルゴリズムによる意思決定の公平性(Fairness)分析や、経済政策の効果測定など、社会現象のメカニズム解明に応用されています。因果グラフは、差別やバイアスに関する仮説を形式化し、検証するための強力なツールとなります。

これらの応用事例に共通するのは、BNが要因間の関係性を視覚化し、データサイエンティストとドメイン専門家の間のコミュニケーションを促進し、データに基づいた意思決定の質を高めることに貢献している点です。

pgmpyの概要

ベイジアンネットワークの理論的な基礎を理解した上で、実際にデータ解析を進めるために不可欠なのが、Pythonライブラリの存在です。pgmpy(Probabilistic Graphical Models in Python)は、確率的グラフィカルモデル(PGM)全般を扱うための、包括的で基礎的なツールキットとして設計されています。

pgmpyは、ベイジアンネットワークの他にも、動的ベイジアンネットワーク(DBN)、マルコフネットワーク、構造方程式モデリング(SEM)といった様々なモデルの構築、学習、推論をサポートしています。このライブラリは、モジュール性と拡張性に焦点を当てた設計思想を持っており、研究者や開発者が基本的な構成要素からモデルを柔軟に構築し、新しいアルゴリズムを容易に拡張・実装できる点が特徴です。

特に、PGMの厳密な理論を学びながら、その実装を深く理解したいユーザーにとってpgmpyは最適な選択肢です。チュートリアルには、Judea Pearl氏の著作に登場する古典的な因果推論の例が豊富に含まれており、教育的な価値も非常に高いと評価されています。

pgmpyの主要な機能とアルゴリズム

pgmpyは、ベイジアンネットワークを用いた因果探索と因果推論のワークフロー全体をサポートする幅広い機能を実装しています。

1. モデルタイプ

pgmpyでは、扱うデータ(確率変数)の性質に応じて、異なるモデルクラスと条件付き確率分布(CPD)が提供されています。

  • 離散ベイジアンネットワーク (DiscreteBayesianNetwork):
    離散的またはカテゴリカルな変数を扱うためのクラスです。CPDのパラメータ化には、確率テーブルとして表現されるTabularCPD(条件付き確率分布表, CPT)を使用します。
  • 線形ガウスベイジアンネットワーク (LinearGaussianBayesianNetwork):
    全ての変数が連続的であり、全てのCPDが線形ガウス分布であるベイジアンネットワークを扱います。このモデルのCPDにはLinearGaussianCPDを使用し、線形回帰の係数と誤差の標準偏差を定義します。

2. 因果探索(構造学習, Structure Learning)

観測データから因果構造であるDAGを推定する機能が充実しています。これは、Causal Discoveryとも呼ばれます。

  • 制約ベースアルゴリズム:
    変数間の条件付き独立性テストの結果を用いてグラフの構造を推定します。代表的なものとして、PCアルゴリズムが実装されており、安定版(stable)や並列版(parallel)などのバリアントも利用可能です。
  • スコアベースアルゴリズム:
    各グラフ構造を評価するためのスコア関数(例:BIC, BDeuなど)を最大化する構造を探索します。HillClimbSearchGES (Greedy Equivalence Search)が主要なアルゴリズムとして提供されています。

3. パラメータ学習 (Parameter Learning)

モデルの構造(DAG)が既知、または構造学習によって得られた後、データを用いてCPDのパラメータを推定します。

4. 因果推論機能

pgmpyのCausalInferenceモジュールは、ベイジアンネットワーク上での介入効果の計算を可能にします。

  • do演算子の実行:
    特定の変数に外部から介入(操作)を行った際の確率分布 \(P(Y|\text{do}(X=x))\) を計算できます。これは、観測された相関 \(P(Y|X=x)\) と因果効果を区別する上で最も重要な概念です。
  • 識別基準のサポート:
    観測データから因果効果を計算するために、バックドア基準(Backdoor Criterion)やフロントドア基準に基づく調整集合の特定と、それを用いた因果効果の識別と計算をサポートしています。この機能は、因果推論の実践において核となる部分です。

pgmpy は、これらの多岐にわたる機能を統合することで、BNを用いた基礎的な学習から高度な因果分析までをカバーする、PythonにおけるPGM研究・開発の標準的な環境を提供しています。

Pythonによる実践

このセクションでは、これまでに学んだベイジアンネットワークの概念を、pgmpy ライブラリを使って実際に動かしながら理解を深めていきます。 データの準備から始まり、因果構造の学習、パラメータの推定、そして最終的な因果推論まで、一連のワークフローをステップバイステップで実践します。

インストール

はじめに、pgmpy ライブラリをインストールします。

$ pip install pgmpy

データセットの準備

まず、分析に使用するデータセットを準備します。ここでは pgmpy に付属しているサンプルモデル ecoli70 を利用します。このモデルからデータをシミュレーションで生成し、それを分析対象の「観測データ」と見なして進めます。

このコードでは、get_example_modelで既存の学習済みモデルをロードし、そのモデルが生成しうるデータをsimulateメソッドで1000件作成しています。これにより、df_continuousというpandas.DataFrameに、後の分析の元となるデータが格納されます。現実の分析では、この部分がCSVファイルやデータベースから観測データを読み込むステップに置き換わります。

# -*- coding: utf-8 -*-

# 必要なライブラリをインポートします
import pandas as pd
from pgmpy.models import DiscreteBayesianNetwork
from pgmpy.estimators import HillClimbSearch, BayesianEstimator, BIC
from pgmpy.inference import VariableElimination
from pgmpy.utils import get_example_model

# --- 1. データセットの準備 ---
# pgmpyにはいくつかのサンプルモデルが用意されています。
# ここでは、大腸菌(E. coli)のタンパク質局在部位に関する'ecoli70'モデルを使用します。
# このモデルは元々、変数間の関係が線形ガウス分布で定義されています。
# 参考: https://www.bnlearn.com/bnrepository/discrete-medium.html#ecoli70
print("--- 1. データセットの準備 ---")
original_model = get_example_model('ecoli70')

# ベイジアンネットワークの学習にはデータが必要です。
# ここでは、既存のモデルからデータを1000件サンプリング(シミュレート)して、学習用データセットを生成します。
# 現実のタスクでは、観測データ(実験結果やアンケートなど)をここに入れることになります。
print("モデルからデータを1000件サンプリング(シミュレート)します...")
df_continuous = original_model.simulate(n=1000, seed=42) # seedを固定して再現性を確保

print("サンプリングした連続値データの最初の5行:")
print(df_continuous.head())
print("-" * 30)

データの離散化

次に、準備したデータをベイジアンネットワークのアルゴリズムで扱える形式に変換します。今回使用する構造学習アルゴリズム(HillClimbSearch)は、離散値を前提としているため、連続値である元のデータをカテゴリデータに変換する「離散化」を行います。

ここでは、pandasqcut関数を用いて、各列(変数)のデータを4つのカテゴリ(0, 1, 2, 3)に変換します。qcutはデータ量を基準に区間を区切る(この場合は四分位数)ため、各カテゴリにほぼ均等な数のデータが割り振られます。

# --- 2. データの離散化 ---
# HillClimbSearchなど、ベイジアンネットワークの多くのアルゴリズムは、
# 「離散値」のデータ(例:「低い」「普通」「高い」)を前提としています。
# 今回のデータは連続値なので、これらをいくつかの区間に区切ってカテゴリに変換します。
#
# ここでは、各変数の値を4つの区間(四分位数)に分け、0, 1, 2, 3 というカテゴリに変換します。
# これにより、例えば「値が小さい方から25%のグループ」がカテゴリ0になります。
print("\n--- 2. データの離散化 ---")
df_discrete = pd.DataFrame()
for col in df_continuous.columns:
    # pandasのqcutは、データを分位数ベースで区切ります。
    # duplicates='drop'は、区切りの値が重複した場合に、区間の数を自動で減らすオプションです。
    df_discrete[col] = pd.qcut(df_continuous[col], q=4, labels=False, duplicates='drop')

print("離散化したデータの最初の5行:")
print(df_discrete.head())
print("(値が 0, 1, 2, 3 などのカテゴリに変換されました)")
print("-" * 30)

因果探索(構造学習)

データの前処理が完了したので、いよいよ因果構造をデータから推定します。このプロセスを「構造学習」または「因果探索」と呼びます。ここでは、数あるアルゴリズムの中から、スコアベース手法の一つである山登り法(Hill-Climb Search) を使用します。

このコードでは、まずHillClimbSearchとスコア関数BICのインスタンスを作成します。そしてhc.estimate()メソッドを実行することで、データに最も適合する(BICスコアが最も高くなる)グラフ構造の探索が開始されます。探索の結果、変数間の依存関係を示すエッジ(矢印)のリストが得られます。これが、データから推定された因果構造の骨格となります。

# --- 3. 因果探索(構造学習) ---
# ベイジアンネットワークの「構造」、つまり変数間の因果関係(矢印の向き)をデータから自動で推定します。
#
# ここでは、山登り法(HillClimbSearch)というアルゴリズムを使用します。
# これは、何もない状態から少しずつ矢印を追加・削除・反転させながら、
# 「スコア」が最も良くなるネットワーク構造を探す、という探索手法です。
#
# スコアリング手法として、BIC(ベイズ情報量規準)を使用します。
# BICは、モデルがデータにどれだけ適合しているか(適合度)と、モデルの複雑さ(辺の数)の
# バランスを取るための指標です。複雑すぎるモデルにはペナルティを与えます。
print("\n--- 3. 因果探索(構造学習) ---")
print("データからネットワーク構造の学習を開始します...")
hc = HillClimbSearch(df_discrete)
scoring_method = BIC(df_discrete)
best_model_structure = hc.estimate(scoring_method=scoring_method)

print("学習したネットワークの構造(変数間の矢印):")
# best_model_structure.edges() は (原因, 結果) のタプルのリストです。
print(sorted(list(best_model_structure.edges())))
print("-" * 30)

モデルの構築とパラメータ学習

因果探索によってグラフの「骨格」が得られたので、次はその骨格に「肉付け」をするパラメータ学習のステップに進みます。具体的には、データを用いて、各変数の条件付き確率分布(CPD)を計算します。これにより、親ノード(原因)が特定の値をとったときに、子ノード(結果)がどの確率で各値をとるかが決まります。

ここでは、ステップ3で得られた構造(best_model_structure.edges())をDiscreteBayesianNetworkクラスに渡してモデルのインスタンスを生成します。その後、model.fit()メソッドに離散化後のデータを渡すことで、パラメータ学習が実行されます。estimatorとしてBayesianEstimatorを指定しており、これにより各CPDがデータに基づいて計算されます。最後にcheck_model()で、構築されたモデルが数学的に整合性が取れているかを確認しています。

# --- 4. モデルの構築とパラメータ学習 ---
# 因果探索で得られた構造をもとに、ベイジアンネットワークモデルの骨格を作成します。
model = DiscreteBayesianNetwork(best_model_structure.edges())

# 次に、各変数が親ノード(原因)の値に応じて、どのような確率で特定の値をとるか(結果)を計算します。
# この確率表を「条件付き確率分布(CPD)」と呼び、これを推定する処理が「パラメータ学習」です。
#
# ここでは、BayesianEstimatorを使用します。これは、データ量が少ない場合でも頑健な推定ができます。
print("\n--- 4. パラメータ学習 ---")
print("ネットワークのパラメータ(条件付き確率)の学習を開始します...")
model.fit(df_discrete, estimator=BayesianEstimator)

# モデルが正しく構築されたか(CPDの確率の合計が1になるかなど)をチェックします
try:
    model.check_model()
    print("モデルの妥当性チェック: 正常です。")
except Exception as e:
    print(f"モデルの妥当性チェック: エラーが発生しました。 {e}")
print("-" * 30)

因果推論

モデルの学習が完了しました。最後に、このモデルを使って因果推論を行います。因果推論とは、ある変数(証拠, Evidence)が特定の値をとったという観測情報が得られたときに、他の変数の確率分布がどのように変化するかを計算することです。ここでは、代表的な推論アルゴリズムである変数消去法(Variable Elimination) を使用します。

このコードブロックでは、まずVariableEliminationクラスに学習済みモデルを渡して推論エンジンを作成します。そしてinference.query()メソッドを用いて、推論を実行します。 例えば、1つ目のクエリP(ygcE | b1191 = 0)は、「b1191という変数がカテゴリ0であると観測された(evidence)場合に、ygcEという変数(variables)が各カテゴリをとる確率はどうなりますか?」という問いに答えるものです。結果として、条件付けられた後のygcEの確率分布が出力されます。

# --- 5. 因果推論 ---
# 学習したモデルを使って、ある事象(証拠)が観測されたときに、
# 他の事象の確率がどう変化するかを計算します。
#
# ここでは、変数消去法(Variable Elimination)というアルゴリズムを使います。
# これは、与えられた証拠に基づいて、知りたい変数の確率分布を正確に計算する代表的な手法です。
print("\n--- 5. 因果推論 ---")
inference = VariableElimination(model)

# 例1: P(ygcE | b1191 = 0) の計算
# 学習された因果関係 ('b1191' -> 'ygcE') に基づいて推論します。
# 「b1191」がカテゴリ `0` であるという情報が得られたとき、「ygcE」の確率分布はどうなるか?
try:
    print("推論クエリ1: P(ygcE | b1191 = 0)")
    query_result1 = inference.query(variables=['ygcE'], evidence={'b1191': 0})
    print("推論結果:")
    print(query_result1)
except Exception as e:
    print(f"推論クエリ1でエラーが発生しました: {e}")

print()

# 例2: P(lpdA | cspG = 1) の計算
# 学習された因果関係 ('cspG' -> 'lpdA') に基づいて推論します。
# 「cspG」がカテゴリ `1` であるとき、「lpdA」の確率分布はどうなるか?
try:
    print("推論クエリ2: P(lpdA | cspG = 1)")
    query_result2 = inference.query(variables=['lpdA'], evidence={'cspG': 1})
    print("推論結果:")
    print(query_result2)
except Exception as e:
    print(f"推論クエリ2でエラーが発生しました: {e}")

print("-" * 30)

おわりに

今回は、ベイジアンネットワークの理論的基盤と、Pythonライブラリ pgmpy を用いた因果探索および因果推論の実践的な側面について解説しました。

ベイジアンネットワークは、単なる予測モデルとしてではなく、データの背後にあるメカニズムを解明するための分析ツールとして活用されています。有向非巡回グラフ(DAG)を用いて因果的な仮定を明示的に表現し、その仮定の下で数学的に正当な推論を導出する理論的基盤を提供します。

この因果関係の理解は、今後のAIシステムに不可欠な要素である透明性(Transparency)、公平性、そして頑健性(Robustness)を高める上で不可欠です。グラフ構造により要因間の関係性が視覚化されるため、モデルが「ホワイトボックス」なAIとして機能し、推論の根拠を説明できるようになります。これは、特に解釈可能性が極めて重要となる実世界の問題、例えば医療診断支援 や品質・歩留まり改善、離職防止施策の立案 など、意思決定の理由を説明する必要がある分野で大きな利点をもたらします。

そして、pgmpyは、このベイジアンネットワークを Python 環境で柔軟に実装できる包括的なツールキットです。因果探索(Structure Learning)やパラメータ学習、そして \(do\) 演算子 を用いた因果推論といった一連のプロセスを統合してサポートし、研究と実務の橋渡しをしています。

More Information