非ガウス性の統計的因果推論

いつか調べるとしていた因果分析ツールですが、統計的因果推論について調べてみることにしました。

1 LiNGAM
2 DirectLiNGAM
3 まとめ
4 参考文献

LiNGAM

LiNGAMは、Linear Non-Gaussian Acyclic Modelの略称です。観測されたデータに対して、何の事前知識もなしに因果構造を発見する方法で、近年では因果探索（Causal Discovery）と呼ばれているようです。[1]

前提条件

データの事前知識は必要ありませんが、いくつかの仮定を満たす必要があります。

観測変数が因果順に並んでいること
- ＝リストの後方の変数が、前方の変数の原因にはならない
- ＝再帰構造にならない
- ＝因果グラフは、Directed Acyclic Graphs（DAG）になる
ある変数は、それより前方の変数の線形関数で表される
- $x_j = \sum_{k(i)<k(j)} b_{ij}x_i + e_i + c_i$
- $e_i$ は、撹乱、誤差、あるいはノイズ
- $c_i$ は、任意定数
- $k(i)$ は、I番目の変数の因果順序の位置
撹乱の特性
- ゼロ分散ではない非ガウス分布に従う連続値乱数
- 互いに独立： $p(e_1,e_2,\cdots)=\prod_{i}p(e_i)$

これらの仮定は、観測されていな交絡因子（cofounders）が無いこと（因果的十分性）を暗黙に仮定しています。

方程式

２番目の線形関数であるという仮定は、次のようなベクトル行列形式で書くことができます。

$\mathrm{\bf x} = \mathrm{\bf Bx}+\mathrm{\bf e}$

そして、この方程式は連立線形方程式になっているので、次のように解くことができます。

$\mathrm{\bf x} = \left(\mathrm{\bf I}-\mathrm{\bf B}\right)^{-1}\mathrm{\bf e}=\mathrm{\bf Ae}$

これは、「解法アルゴリズムの選び方」で説明した通り、解を線型結合（線形関数）に無理やり置き換えたことで、因果関係を連立線形方程式問題に置き換えたと見ることができます。

ここで、行列Bは、因果順の仮定により、対角成分ゼロの下三角行列になります。そのため、行列I–Bは、対角成分１の下三角行列にになります。これは、LU分解の行列Lと同じ形です。そして、この形の行列の場合、ガウスの消去法で必ず解を求めることができます。逆行列Aも、うまく順序を入れ替えると、対角成分がある下三角行列になります。

解法

さて、問題 $\mathrm{\bf x} = \mathrm{\bf Ae}$ は、観測データxが既知で、混合行列Aと撹乱eが未知の場合、線形の独立成分分析（Independent Component Analysis : ICA）の問題に等しくなります。この問題は、十分に多い観測データxがあると、混合行列Aを推定できることが示されています。ここで、独立成分とは、誤差eのことです。

ただし、ICAは因果関係にとって不都合な点が２つあります。それは次のような点です。

ICAでは独立成分の順序は任意だが、因果関係では因果順に並んでいなければならない（仮定１）
ICAでは行列Aを適切にスケーリングするが、因果関係では対角成分が１に固定されていなければならない

そのため、因果関係を推定するには、ICAを計算する→規格化する→因果順に置換する、という手順が必要です。

アルゴリズム

以上により、仮定１〜３の下では、次のような因果発見アルゴリズムを構築できます。

LiNGAMアルゴリズム

データ行列X（m×n）を、ICAによりX=ASと分解し、W=A^{-1}とする。（m:成分数、n:サンプル数、m<<n、S:撹乱行列（m×n）、A:混合行列（m×m））
Wの成分を並べ替え、非ゼロ対角成分を持つW’を作成する。（ゼロ割防止）
W’の各行をその行の対角成分で割り、対角成分を１にする。（規格化）
B’=I–W’を計算する。
B=PB’P^Tが下三角行列となるような置換行列Pを見つけるため、 $f\left(\mathrm{\bf P}\right)=\sum_{i\ge j}B_{ij}^2$ を最小化する。（因果順）

ステップ５の最適化問題の目的関数f(P)は、下三角行列Bの対角成分および上三角成分はゼロにならないといけないことに由来します。各成分を二乗することで正負成分の打ち消しを防ぎ、二乗和f(P)=0の時の置換行列Pを探索する最適化問題にしています。（参考記事）

アルゴリズムのより詳細な説明は、参考文献１をご覧ください。

DirectLiNGAM

現在では、上記のICAを使ったLiNGAMではなく、DirectLiNGAMの方がデフォルトになっているようです。

ICAを使ったLiNGAMには、２つの問題点がありました。[2]

ICAアルゴリズムで行われる最適化が、パラメータの設定のよっては収束しない。または、収束するパラメータを見つけることが難しい。
置換アルゴリズムがスケール不変ではなく、スケールによっておかしな因果順が推定されてしまう。しかし、本来、因果順はスケールに依存しないはずである。

これらを解決するために、DirectLiNGAMではICAを使わずに因果探索を行います。

アプローチ

DirectLiNGAMの計算アプローチは、下図のように、補題１・補題２・系１の３つの矢印をたどることで計算されます。

LiNGAMでは入力データ $x_i$ について因果探索を行いましたが、DirectLiNGAMでは入力データの残差 $r_i^{(j)}$ について因果探索を行う点が大きく異なります。

DirectLiNGAMのアプローチでは、ICAは、入力データから残差への変換に置き換えられます。

理論

補題１は、次のようなものです。

補題１
インプットデータ $\mathrm{\bf x}$ が、厳密にLiNGAM条件（３つの前提とデータの有限性）を満たすとする。 $x_i$ が $x_j$ について回帰された残差を $r_i^{(j)}=x_i-\displaystyle\frac{\mathrm{cov}(x_i,x_j)}{\mathrm{var}(x_j)}x_j,\ \left(i\neq j\right)$ と書くことにする。このとき、 $x_j$ が、全ての残差 $r_i^{(j)} (i = 1,2,\cdots,j-1,j+1,\cdots)$ と独立であれば、 $x_j$ は外生変数である。
Shimizu, S., et. al (2011)

外生変数とは、因果の連鎖の出発点になる変数のことで、出力の矢印しかない変数のことです。言い換えると、入力データの中には説明変数となるデータが存在しない変数のことで、入力されたどのデータにも従属しない変数（＝独立変数）となります。

補題１は、あるデータ $x_j$ が、それ以外のデータとの回帰残差 $r_i^{(j)}$ が全て独立なら、独立変数であることを示しています。ただし、重要なのは逆の場合で、データ $x_j$ と残差 $r_i^{(j)}$ が非独立の場合、 $x_i\rightarrow x_j$ という従属関係があることになります。

この非独立性が成立するには、撹乱 $e_i$ の非ガウス性（ガウス分布にならないこと）が必要です。これは、Darmois-Skitovitchの定理に依拠しています。

Darmois-Skitovitchの定理
ランダムな独立変数 $s_i, (i=1,2,\cdots,q)$ の線型結合によって、２つのランダムな変数 $y_1, y_2$ を、 $y_1=\displaystyle\sum_{i=1}^q\alpha_i s_i,\quad y_2=\displaystyle\sum_{i=1}^q\beta_i s_i$ と定義する。このとき、 $y_1,\ y_2$ が独立ならば、 $\alpha_j\beta_j\neq 0$ となる全ての変数 $s_j$ はガウス分布となる。
Shimizu, S., et. al (2011)

この定理は、逆に解釈すると、変数 $y_1,\ y_2$ が非独立ならば変数 $s_j$ は非ガウス分布になることを示しています。ここで、 $y_1,\ y_2,\ s_j$ が、 $x_j,\ r_i^{(j)},\ e_l$ となれば、非独立性のために撹乱の非ガウス性が必要なことを示すことができます。

証明方針は、外生変数は $x_i=e_i$ と撹乱そのものになることを使って、（１）全ての従属関係を代入すれば任意の変数は $x_j=f(\mathrm{\bf e})$ と撹乱のみで表示できること、（２）このように外生変数のみで表示された変数を代入すれば残差も $r_i^{(j)}=g(\mathrm{\bf e})$ と撹乱のみで表示できること、（３）この二つにDarmois-Skitovitchの定理を適用することで、非独立なら非ガウスを示すことができます。

詳しい証明は、参考文献２をご覧ください。

補題２とその系１は、次のようなものです。

補題２
インプットデータ $\mathrm{\bf x}$ が、厳密にLiNGAM条件（３つの前提とデータの有限性）を満たすとする。さらに、 $x_j$ は外生変数だとする。ここで、 $x_j, (i\neq j)$ で回帰された(p-1)個の $x_i$ の残差 $r_i^{(j)}$ を集めた(p-1)次元の残差ベクトルを $\mathrm{\bf r}^{(j)}$ と書くことにする。このとき、残差ベクトル $\mathrm{\bf r}^{(j)}$ は、 $\mathrm{\bf r}^{(j)}=\mathrm{\bf B}^{(j)}\mathrm{\bf r}^{(j)}+\mathrm{\bf e}^{(j)}$ を満たす。ここで、 $\mathrm{\bf B}^{(j)}$ は行と列の同時並べ替えで厳密に下三角になるように並べ替えられる行列で、 $\mathrm{\bf e}^{(j)}$ の要素は非ガウス分布かつ互いに独立である。
Shimizu, S., et. al (2011)

系１
インプットデータ $\mathrm{\bf x}$ が、厳密にLiNGAM条件（３つの前提とデータの有限性）を満たすとする。さらに、 $x_j$ は外生変数だとする。残差 $r_i^{(j)}$ の因果順序を $k_{r^{(j)}}(i)$ と書き、外生変数 $x_i$ の因果順序を $k(i)$ と書くことにする。このとき、残差 $r_i^{(j)}$ の因果順序は、変数 $x_i$ の因果順序は同じになる。すなわち、 $k_{r^{(j)}}(l)< k_{r^{(j)}}(m) \Leftrightarrow k(l)< k(m)$ となる。
Shimizu, S., et. al (2011)

直観的には、データxがLiNGAM型の方程式x=Bx+eを満たす場合、データの線型結合である回帰残差rもr=Br+eを満たすことは想像できると思います。

ですが、補題２と系１の証明はいまいち理解できませんでした。もし、理解できる方がいたら教えていただけると助かります。

アルゴリズム

これらを踏まえると、次のようなアルゴリズムを組むことができます。

DirectLiNGAMアルゴリズム

初期化
- p次元乱数ベクトル $\mathrm{\bf x}$
- 添字集合 $U=\{1,2,\cdots,p\}$
- p×n次元乱数行列 $\mathrm{\bf X}$
- 空のカーネル添字集合 $K=\{\phi\}$
- 外生変数の固定添字 $m:=1$
For j=1…p-1:
1. Kに含まれずUに含まれるの添字iについて、 $x_j$ に対して $x_i$ を最小二乗法で回帰計算し、残差ベクトル $\mathrm{\bf r}^{(j)}$ を作る。
2. 同様に、データ行列 $\mathrm{\bf X}$ から残差行列 $\mathrm{\bf R}^{(j)}$ を作成する。
3. 残差ベクトル $\mathrm{\bf r}^{(j)}$ に最も独立したデータ $x_m$ を探索する。
4. 独立データの添字mを、Kに追加する。
5. $\mathrm{\bf x}=\mathrm{\bf r}^{(m)}$ 、 $\mathrm{\bf X}=\mathrm{\bf R}^{(m)}$
残ったデータの添字をKに追加する。（Kは添字が因果順序に並んでいる）
Kの順序で下三角行列Bを構成し、元のデータ $\mathrm{\bf x}$ や $\mathrm{\bf X}$ を最小二乗法や最尤法で回帰して、行列Bの要素を計算する。