はじめに
変数間の因果関係を推測する学問として、統計的因果探索というものがあります。その中でLinGAMというモデルについて説明します。
まず始めにLinGAMとその推測方法について説明して、最後にpythonで実際にLinGAMを使ってみます。
使えればいいよと言う方は最後だけご覧ください。
統計的因果探索とは？
はじめに、簡単に統計的因果探索について説明します。
統計的因果探索は一言で言うとデータから因果関係を推測することです。
因果関係とは、原因と結果のことで例えば、
勉強時間とテストの成績には相関があることがわかったとします。この場合どちらが原因でどちらが結果になるでしょうか？
もちろん勉強時間が原因で成績が結果になります。
このような原因と結果の因果関係をデータから推測することを統計的因果探索と呼び、今回のブログでLinGAMと言うモデルを用いた統計的因果探索手法について説明します。
LinGAM
LinGAMとはLinear Non Gaussian Modelのことで線形モデルで非ガウス分布を仮定したモデルのことを言い、以下のように定式化されます。
x = Ax + \epsilon
x = (x_1~x_2...x_n)^T
\epsilon = (\epsilon_1~\epsilon_2...\epsilon_n)^T
ここで \epsilon_i はそれぞれ独立に非正規分布に従い、Aの対角成分は0かつ上三角行列にできる。
対角成分が0であるのはx => xと言う因果関係をなくすためです。
また上三角行列の仮定は、因果関係が非循環であることを意味しています。
なぜ非正規分布？
なぜ非正規分布を仮定しているのでしょうか？ \epsilon が正規分布に従う場合と、一様分布に従う場合を比較してみましょう。
x = y + \epsilon と y = x + \epsilon について比較します
import numpy as np
import matplotlib.pyplot as pltまずは正規分布の場合を考えます。
注) 共分散が等しくなるようにしてあります。
# 正規分布から発生してるモデルを二つ考える
# X1 => Y1のモデル
X1 = np.random.randn(5000)
Y1 = X1*0.8 + np.random.normal(0,0.36,5000)
# Y2 => X2のモデル
Y2 = np.random.randn(5000)
X2 = Y2*0.8 + np.random.normal(0,0.36,5000)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.scatter(X1,Y1)
plt.title('X1 => Y1')
plt.subplot(1,2,2)
plt.scatter(X2,Y2)
plt.title('Y2 => X2')
plt.show()
上記からもわかるように、二つの分布は等しくなっており、データから因果関係を区別することはできなくなっています。
次に一様分布の場合を考えます。こちらも同様に共分散が等しくなるように設定してあります。
# 一様分布から発生してるモデルを二つ考える
# X1 => Y1のモデル
X1 = np.random.uniform(-12**(0.5)/2,12**(0.5)/2,5000)
Y1 = X1*0.8 + np.random.uniform(-0.6*12**(0.5)/2,0.6*12**(0.5)/2,5000)
# Y2 => X2のモデル
Y2 = np.random.uniform(-12**(0.5)/2,12**(0.5)/2,5000)
X2 = Y2*0.8 + np.random.uniform(-0.6*12**(0.5)/2,0.6*12**(0.5)/2,5000)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.scatter(X1,Y1)
plt.title('X1 => Y1')
plt.subplot(1,2,2)
plt.scatter(X2,Y2)
plt.title('Y2 => X2')
plt.show()
共分散は等しいですが、上記の二つの分布は異なっていることがわかります。
つまり、データから因果関係を推測することが可能であることがわかります。また、一様分布以外でも非正規分布では上記同様のことが成り立ちます。
推測方法
非正規分布を仮定した場合、共分散が等しい(等しい相関を持つ)場合でも因果関係を推測できることが確認できました。
では、実際にどのように推測するのかについて説明します。
ここでは前回ブログに書いた独立成分分析を利用します。
ここで改めてLinGAMモデルを確認します。
x = Ax + \epsilon
x = (x_1~x_2...x_n)^T
\epsilon = (\epsilon_1~\epsilon_2...\epsilon_n)^T
次のように式変形をします。
x = (I - A)^{-1}\epsilon
\epsilon_iがそれぞれ独立であることから独立成分分析のモデルに等しくなっていることがわかります。
つまり、独立成分分析により、 (A-I)^{-1}の尺度と行の交換を除きPD(I-A)^{-1}が推測できることになります。
Pは行の交換を行う行列、Dは尺度を変更する対角成分のみ値を持つ行列であり、
独立成分分析ではP、Dは特定できません。
尺度と行の交換
独立成分分析ではP、Dは特定できませんが、LinGAMでは特定することができます。
なぜなら、Aの対角成分が0であり上三角行列にできるという仮定があるからです。この時I - Aの対角成分は1になります。
(PD(I-A)^{-1})^{-1} = (I-A)D^{-1}P^{-1}
となるのでPD(I-A)を独立成分分析で見つけた後Aが仮定を満たすようにP、Dを決定することでAを決定することができます。
このようにしてLinGAMでは因果関係を推測することができます
pythonで実際に因果関係を推測する。
先ほど例で一様分布から発生させたサンプルで試してみます。
! pip install lingam
import lingamdata1 = np.concatenate((X1.reshape(-1,1),Y1.reshape(-1,1)),1)
data2 = np.concatenate((X2.reshape(-1,1),Y2.reshape(-1,1)),1)
label1 = ['X1','Y1']
label2 = ['X2','Y2']
lg = lingam.DirectLiNGAM()
lg.fit(data1)
causal_order = lg.causal_order_
B1 = lg.adjacency_matrix_
print('因果順序は {} => {} '.format(label1[causal_order[0]],label1[causal_order[1]]))
print('回帰係数は')
print(B1)
print('')
print('')
lg = lingam.DirectLiNGAM()
lg.fit(data2)
causal_order = lg.causal_order_
B2 = lg.adjacency_matrix_
print('因果順序は {} => {} '.format(label2[causal_order[0]],label2[causal_order[1]]))
print('回帰係数は')
print(B2)因果順序は X1 => Y1
回帰係数は
[[0. 0. ]
[0.80440109 0. ]]
因果順序は Y2 => X2
回帰係数は
[[0. 0.80257153]
[0. 0. ]]適切に推測が行えていることがわかります。
では次に3変数の場合で検証してみます。
x1 = np.random.uniform(-1,1,1000)
x2 = x1*0.5 + np.random.uniform(-1,1,1000)
x3 = x1*-0.5 + x2 + np.random.uniform(-1,1,1000)
data = np.concatenate((x1.reshape(-1,1),x2.reshape(-1,1),x3.reshape(-1,1)),1)
lg = lingam.DirectLiNGAM()
lg.fit(data)
causal_order = lg.causal_order_
B = lg.adjacency_matrix_
print('回帰係数は')
print(B)
回帰係数は
[[ 0. 0. 0. ]
[ 0.49574869 0. 0. ]
[-0.45706713 0.98177538 0. ]]3変数の場合でも適切に回帰係数が計算できていることがわかります。
参考文献
清水 昌平 , 統計的因果探索(機械学習プロフェッショナルシリーズ), 講談社,2017
