統計的因果探索の1つであるLinear non-Gaussian Acyclic Model(LiNGAM)をpython3で書いてみました。githubに置いてます→ICA-LiNGAM code
ちょっと説明
どんなモデル?
LiNGAMでは、データの生成過程を因果と捉える手法で、仮定に沿っているデータであれば理論的には因果関係を一意に推定可能です！(線形関係)
要は、下のような状況で生成されているデータがあったとしたら、
このデータを与えるだけで、"xからyに強さ2の因果効果がある"と分かる手法です。
x = e1
y = 2x + e2
すごさ
何がすごいかというと、重回帰分析とかだとyをxで説明するとか、xをyで説明するとかの方向を自分で決めなければいけなかったのに、LiNGAMではデータから確実に決められるところが画期的です。
他にも仮定の自然さみたいなのが挙げられます。
その仮定は、実データを考えてみると結構自然なもので、
1. データが連続変数
2. ノイズ項eが、非ガウス分布(正規分布でない)、
3. ノイズ項eが、互いに独立
であるということです。ちなみにノイズ項は外生変数とか言われます。
LiNGAMは、非ガウス性に注目してるおかげで、連続変数同士の因果関係を見つけることができる手法です。
実装
内部で独立成分分析(ICA)という信号処理を使っていて、sklearnにFastICAがあったので、割と簡単にできました。
ただ、非ガウス性の尺度として、ネゲントロピー(近似)を扱ってるものしか発見できなかったので、今回は尖度ベースのICAを自分で実装してみました。
コードは全てココにあります。
実行例
Step1
まず、lingam.pyを実行する。jupyterだと%run lingam.pyとかでできる。
Step2
サンプルデータ生成する。この例では、x -> y の向きに2の強さの因果がある。
一様分布で非ガウスなノイズを生成。
import numpy as np
import pandas as pd
size = 10000
np.random.seed(2017)
x = np.random.uniform(size=size)
np.random.seed(1028)
y = 2*x + np.random.uniform(size=size)
X = pd.DataFrame(np.asarray([x,y]).T,columns=["x","y"])
Step3
実行する。fitでできる。
lingam = LiNGAM()
lingam.fit(X)
lingam = LiNGAM()
lingam.fit(X,use_sklearn=True)
で、実行結果。この例ではどっちも同じだった。
ほぼ正しく因果が推定できてる!!
x ---|1.991|---> y
array([[ 0. , 0. ],
[ 1.99149033, 0. ]])
おわり
3以上変数でもできます。ただ、あんまり多いとすごく時間かかるかも...
尖度使ったICAは、異常値に弱いので、実データならばuse_sklearn=Trueの方がいいと思います笑
質問、ミスなどありましたら、教えてください！
他の因果探索手法も実装してみたい！