統計的因果探索：ノート4章

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。

統計的因果探索 (機械学習プロフェッショナルシリーズ)

作者:清水昌平
発売日: 2017/05/25
メディア: 単行本（ソフトカバー）

著者のスライド

まとめ

4章：実際に LiNGAM で変数間の因果関係を推定する方法は（この本で主に取り上げられているのは）2つあって、
- 観測変数間の関係を全体的に求めてそれと整合的な因果グラフを特定する派：
  まず $x = (I-B)^{-1}e$ の $I-B$ を求める（ICA の手法でスケーリングと行の順序を除いて特定できる）
  → 「対角成分が1」「下三角」を満たすようにスケーリングと行の順序を特定する
  → ここまでで観測変数ベクトルの成分が因果的順序になっているので後は自分より親で回帰
- 2観測変数間の関係からどの観測変数が祖先が順に特定していく派：
  因果が正しい向きに回帰すれば誤差と説明変数が独立になることを利用し親のない観測変数を特定する
  → 全体からその変数の寄与を除いて繰り返す

前回：ノート1章、2章、3章

f:id:cookie-box:20200101101603p:plain:w60

3章までで構造的因果モデルに線形関数・非ガウスノイズを仮定すれば観測変数の分布から因果グラフを識別できることはわかりましたが、その具体的な手法はまだ示されていませんね。これでは、ありうるすべての因果グラフにありうるすべての誤差分布を試さなければならないのではないでしょうか。早速4章を読んで LiNGAM なる手法の詳細を学びましょ…独立成分分析（ICA）？最初から LiNGAM を教えてくれるのではないのですか？仕方ないので読んでいくよりないですね…以下のようなモデルを考えるとのことです。

$\begin{split} x_1 &= a_{11} s_1 + a_{12} s_2 \\ x_2 &= a_{21} s_1 + a_{22} s_2 \end{split}$ 各変数の意味は例えば以下です。

$x_1$ ：マイク1が集音した波形（のあるフレームの値）
$x_2$ ：マイク2が集音した波形（のあるフレームの値）
$s_1$ ：人物1の声の波形（のあるフレームの値）
$s_2$ ：人物2の声の波形（のあるフレームの値）

それで、ICA の目的はデータ行列 $X$ （ここでは、1行目にマイク1が集音した波形を、2行目にマイク2が集音した波形を $n$ フレーム並べた $2 \times n$ 行列ですね）から $a_{11}, \, a_{12}, \, a_{21}, \, a_{22}$ を推定して $s_1, \, s_2$ を復元すること…なんと、ということはもしや、ユニット曲から個々のアイドルのボーカルを抽出できるのですか？

f:id:cookie-box:20200101101614p:plain:w60

違うんじゃないかな。88ページに ICA の仮定として $s_1, \, s_2$ が独立であるとあるよ。同じ曲を歌う2人の歌声は独立にならないでしょ。その曲の楽譜という共通原因が存在するよ。もしかしたらそれを独立になるようにプレ処理するフレームワークとかがあって歌声にも適用できるのかもしれないけど知らない。…本題に戻ると、ICA の目的は統計的因果探索の目的ととても似ているね。 $X$ から構造方程式の係数行列を推定するところまでは。でもそこから先が、ICA は「未観測変数を復元すること」、統計的因果探索は「因果グラフを復元し平均因果効果を推定すること」とちょっとずれるね。あと、ICA では観測変数が他の観測変数の原因になることがないのかな。マイク1の波形がマイク2の波形に因果効果をもたらすことはない。

$s_1, \, s_2$ は独立なのですか。2つのマイクの側で2人の人間が互いを気にせずにしゃべっているとはどういう状況なんです。まあいいです。 $s_1, \, s_2$ へのもう1つの仮定は、これらが非ガウス分布にしたがうということですね。それで、 $x = As$ の $A$ は正則と仮定するんですね。もし上の例の 2×2 行列 $A$ のランクが 2 ではなくて 1 だったら、例えば $A = \left( \begin{array}{cc} 1 & 2 \\ 2 & 4 \end{array} \right)$ だったら、マイク1もマイク2もどちらも人物1より人物2の声の方を2倍強く集音する全く同じ役割のマイクになっており、マイクを2つ用意した意味がありません。予算の無駄すぎます。こんなミスが起きていないと信じられる状況であれば「 $A$ はわからないけど正則ではあるはずだ」というのは無理な仮定ではないのでしょうかね。まあ認めましょう。この仮定をすれば、「 $A$ を一意に推定することはできませんが、それと列の順序と尺度だけが異なる可能性のある行列なら推定できるということです（89ページ）」ということです。はい？尺度とは？？

その後を読めばわかる。 $A$ と列の順序と尺度だけが異なる可能性のある行列 $A_{\rm ICA}$ は、対角成分がゼロでない対角行列 $D$ と置換行列 $P$ を用いて $A_{\rm ICA} = ADP$ とかけるということだ。これは $A$ の $j$ 列を $d_{jj}$ 倍（ $d_{jj} \neq 0$ ）してから列をシャッフルする操作だね。ICA で特定できる $A_{\rm ICA}$ は本当にほしい $A$ に右から $D$ と $P$ がかかってしまっているということだ。 $A$ の各列のベクトルの向きだけが特定できると言い換えてもいいだろう（ $d_{jj}$ が負かもしれないからベクトルの向きを特定できるというのもちょっと語弊があって、ベクトルがのる直線を特定できるといった方がいいのかな）。なぜそんな中途半端にしか特定できないのかは理解しやすい。まず、人物1の声を第1成分にして人物2の声を第2成分にする必然性がない。逆にしても構造方程式はかける。そして、人物1の声を2倍にして人物1の声にかかる係数をすべて 1/2 倍しても構造方程式はかける。ここも必然性がない。…この状況で推定できる可能性があるのは、1人目（人物1か人物2かはわからない）の声を定数倍したものと、2人目の声を定数倍したものだね。なるほど「独立成分分析」だ。

92ページには、もしガウス分布だったら $A_{\rm Gaussian} = ADPQ$ までしか特定できないとありますね。 $Q$ は直交行列です。t分布の話のときに出てきました。多変量標準正規分布にしたがう確率ベクトルを $Q$ で変換しても多変量標準正規分布にしたがう確率ベクトルです。t分布の話のときはこれを利用してもとの確率ベクトルの各成分の平均を絞り出しました。

各成分から自由にちょっとずつ絞り出せちゃったら困っちゃうんだよね。そんな変換をしたベクトルはもう、人物1の声も人物2の声も混ざり合っていて分離できてない。

$x=($ 正則な係数行列 $)($ 各成分が独立な確率分布にしたがうベクトル $)$ という形に分解してくださいといわれたとき、 $x=(A)(s)$ と $x=(ADPQ)(Q^\top P^{-1} D^{-1} s)$ を区別できないということですね、確率分布がガウス分布であると仮定すると。…まあそれで、 $A$ を推定するには $y = Wx$ の成分どうしが独立になるような $W$ を見つけて逆行列をとるんですね。具体的には、 $y$ の相互情報量を最小にする…相互情報量は、その多変量確率分布と、各成分どうしが独立であると考えた場合の確率分布のカルバック・ライブラー情報量ですよね。相互情報量=0 を達成できたらそのときはもう各成分どうしは独立です。…各成分どうしが独立の分布とのカルバック・ライブラー情報量とはどうも既視感がありますね…平均場近似です。しかし、あのときはカルバック・ライブラー情報量の向きが逆でしたね…。

平均場近似の文脈では「各成分どうしが独立の分布」は「手段」だったけど、今回は「目的」だしね…。

平均場近似 ― データの生成分布を近似的に推定するために、独立な成分（のブロック）が2つ以上あるモデルを仮定し、交互に更新する。
独立成分分析 ― データの独立成分を復元する。

まあそれで、データから相互情報量を推定する方法は、いまは観測の次元数 $p$ が復元したい独立成分の数 $q$ に等しいとして、データ行列 $X$ に $W \in \mathbb{R}^{q \times q}$ をかけて独立成分に分解しようとしているとして、
f:id:cookie-box:20200223181619p:plain:w630

f:id:cookie-box:20200223181619p:plain:w630

以下を最小にする $W$ を推定するらしいんだけど、とするとこの $p(\cdot)$ が何なのかな…その $W$ のときの $y$ の分布を推定する？具体的には不動点法とかのアルゴリズムがあるらしい…。
$\hat{I}(y) = \displaystyle \left\{ \sum_{j=1}^q \frac{1}{n} \sum_{m=1}^n - \log p(y_j^{(m)}) \right\} - \frac{1}{n} \sum_{m=1}^{n} - \log p(y^{(m)})$

ともかく LiNGAM に進みましょう。構造方程式は $x_i = \sum_{j \neq i} b_{ij} x_j + e_i \, (i = 1, \cdots p)$ ですね。行列表現すれば $x = Bx + e$ です。96ページには非巡回性の話がありますね。 $k(i)$ というのは「 $x_i$ が $x_j$ の原因になるならば必ず $k(i) < k(j)$ となるような各成分への番号の振り方」ですよね。添え字 $i$ の成分が添え字 $k(i)$ になるような並べ替えを行えば $B$ が狭義下三角になると。 $k(i)$ が定義できれば $B$ が狭義下三角になるのはわかりますが、 $k(i)$ が常に定義できるんでしょうか…98ページの下から5行目に「必ずあります」とありました。

証明の仕方はわからないけど、閉路がないグラフなら、「上流」から値を流していったときに2回以上値が更新されるノードはないわけで、すべてのノードが「移項」をしなくてもすんなり値が求まるわけで、それってガウスの掃き出し法でいう「前進消去」が既にできている状態だと思う。だから並べ替えさえちゃんとすれば $B$ が狭義下三角（ $I-B$ が対角成分が1の下三角）になるって気はするかな…。

ともあれ、 $B$ は狭義下三角としてよいのですね。それで、 $x = Bx + e$ $\; \Leftrightarrow \;$ $x = (I-B)^{-1}e \equiv Ae$ と変形して…結局独立成分分析じゃないですか！いやでも待ってください、 $I-B$ って必ず正則になりますか？

必ず正則になるよ。 ${\rm det}(I-B) = 1 \neq 0$ だからね。

上三角行列と下三角行列の性質 | 高校数学の美しい物語

う、確かに。いやだとしても、ICA では $A$ を $ADP$ までしか特定できませんでしたよね。となると $(I-B)^{-1}$ も $(I-B)^{-1}DP$ までしか特定できず、結局 $I-B$ は $P'D'(I-B)$ までしか特定できず、因果グラフが定まらないじゃないですか。

違うかな。101～105ページを読めばわかる。つまり、ICA では $A$ は正則であるという情報しかなかったけど、いま LiNGAM では $I-B$ は「対角成分が1の下三角行列である」までわかっている。であれば、求まった $P'D'(I-B)$ に、下三角にするような行の入れ替え $P'^\top$ と、対角成分を1にするようなスケーリング $D'^{-1}$ を順に左からかけてやればいい。 $I-B$ は特定できる。理論的にはね。

なんと、ぴったり復元できますね。人生で一番感動しました…。

人生で！？それで、4章の続きは具体的な推定アルゴリズムだね。1つ目は ICA の解法を活用する方法で、手順は以下かな。

ICA の何らかの推定アルゴリズムで $\hat{W}_{\rm ICA}$ を求める。
次に、対角成分をゼロでなくすため、対角成分の絶対値を大きくする置換 $\displaystyle \hat{\tilde{P}} = \underset{\tilde{P}}{\rm min} \frac{1}{ | (\tilde{P}\hat{W}_{\rm ICA})_{ii}|}$ を求めて、対角成分が1になるようにスケーリングして $\hat{B}$ を得る。
次に、狭義下三角にするため、 $\ddot{P} x = (\ddot{P} \hat{B} \ddot{P}^\top) \ddot{P} x + \ddot{P} e$ より、 $\ddot{P} \hat{B} \ddot{P}^\top$ が下三角行列に近くなる置換 $\displaystyle \hat{\ddot{P}} = \underset{\ddot{P}}{\rm min} \sum_{i \leqq j} (\ddot{P} \hat{B} \ddot{P}^\top)_{ij}^2$ を求める（ $p$ が大きくなると全探索は困難になるので、 $\hat{B}$ の絶対値が小さい方から $p(P+1)/2$ 個の成分をゼロにしたら狭義下三角にできるかを調べるなどで回避する）。
ここまでで変数間の因果的順序が求まるので、各変数を親変数候補（因果的順序が自分より先である変数の集合）で回帰し、 $B$ を推定する。ただし、普通に回帰したのでは冗長な有向辺が出るので、一致推定量 $\hat{b}_{ij}$ を求めた後適応型 Lasso $\| x_i - \sum_{k(j) < k(i)} b_{ij} x_j \|^2 + \lambda \sum_{k(j) < k(i)} |b_{ij}| / |\hat{b}_{ij}|^\gamma$ する。

手順の3つ目の狭義下三角にするのは以下のようなイメージだね。
f:id:cookie-box:20200224095402p:plain:w320

まず対角成分を確定させるんですね。しかし、対角成分たちとして絶対値の逆数の和が最小な組み合わせを選べば大丈夫なんでしょうか。

どうなんだろ。詳しくかいてないからわかんないな。もう1つの推定アルゴリズムは LiNGAM 専用アプローチだね。

正しい因果的順序で回帰すれば説明変数と残差が独立になるが正しくない順序だと従属になる？なぜこのような差が…慢心、環境の違い…。ときに副部長、説明変数と被説明変数を誤ってしまうと本当に従属になるのですか？

なんで？

以下の（１）（２）が成り立っているかどうかを確認してみますよ。

（１） $x_2$ を $x_1$ で回帰した残差 $\displaystyle r_2^{(1)} = x_2 - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_1)} x_1$ が $x_1$ と独立である。
（２） $x_1$ を $x_2$ で回帰した残差 $\displaystyle r_1^{(2)} = x_1 - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} x_2$ が $x_2$ と独立である。

まず準備として、以下のようになっていますよね。
$\mathbb{E}( x_1 ) = \mathbb{E}( e_1 )$
$\mathbb{E}( x_2 ) = b_{21} \mathbb{E}( e_1 ) + \mathbb{E}( e_2 )$
$\displaystyle \mathbb{E}( r_2^{(1)} ) = b_{21} \mathbb{E}( e_1 ) + \mathbb{E}( e_2 ) - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_1)} \mathbb{E}( e_1 ) = \left( b_{21} - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_1)} \right) \mathbb{E}( e_1 ) + \mathbb{E}( e_2 ) = \mathbb{E}( e_2 )$
$\displaystyle \mathbb{E}( r_1^{(2)} ) = \mathbb{E}( e_1 ) - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \bigl( b_{21} \mathbb{E}( e_1 ) + \mathbb{E}( e_2 ) \bigr) = \left( 1 - b_{21} \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \right) \mathbb{E}( e_1 ) - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \mathbb{E}( e_2 )$

だから（１）については以下のようになって独立になっています。
$\displaystyle \mathbb{E}( r_2^{(1)} x_1 ) = \left( b_{21} - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_1)} \right) \mathbb{E}( e_1^2 ) + \mathbb{E}( e_1 ) \mathbb{E}( e_2 ) = \mathbb{E}( e_1 ) \mathbb{E}( e_2 )$
$\displaystyle \mathbb{E}( r_2^{(1)} ) \mathbb{E}( x_1 ) = \left( b_{21} - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_1)} \right) \mathbb{E}( e_1 )^2 + \mathbb{E}( e_1 ) \mathbb{E}( e_2 ) = \mathbb{E}( e_1 ) \mathbb{E}( e_2 )$
$\Rightarrow \, \mathbb{E}( r_2^{(1)} x_1 ) - \mathbb{E}( r_2^{(1)} ) \mathbb{E}( x_1 ) = 0$

他方、（２）については以下のようになって、やはり独立になっています。
$\displaystyle \mathbb{E}( r_1^{(2)} x_2 ) = b_{21} \left( 1 - b_{21} \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \right) \mathbb{E}( e_1^2 ) + \left( 1 - 2 b_{21} \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \right) \mathbb{E}( e_1 ) \mathbb{E}( e_2 ) - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \mathbb{E}( e_2^2 )$
$\displaystyle \mathbb{E}( r_1^{(2)}) \mathbb{E}( x_2 ) = b_{21} \left( 1 - b_{21} \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \right) \mathbb{E}( e_1 )^2 + \left( 1 - 2 b_{21} \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \right) \mathbb{E}( e_1 ) \mathbb{E}( e_2 ) - \frac{{\rm cov}(x_1, x_2)}{{\rm var}(x_2)} \mathbb{E}( e_2 )^2$
$\Rightarrow \, \displaystyle \mathbb{E}( r_1^{(2)} x_2 ) - \mathbb{E}( r_1^{(2)} ) \mathbb{E}( x_2 ) = 0$

より一般に、因果関係のわからない $x_b$ を $x_a$ で回帰するとき、
$\displaystyle \mathbb{E}( r_b^{(a)} ) = \mathbb{E}(x_b) - \frac{{\rm cov}(x_a, x_b)}{{\rm var}(x_a)} \mathbb{E}(x_a)$
$\displaystyle \mathbb{E}( r_b^{(a)} x_a ) = \mathbb{E}(x_a x_b) - \frac{{\rm cov}(x_a, x_b)}{{\rm var}(x_a)} \mathbb{E}(x_a^2)$
$\displaystyle \mathbb{E}( r_b^{(a)}) \mathbb{E}( x_a ) = \mathbb{E}(x_a) \mathbb{E}(x_b) - \frac{{\rm cov}(x_a, x_b)}{{\rm var}(x_a)} \mathbb{E}(x_a)^2$
$\displaystyle \mathbb{E}( r_b^{(a)} x_a ) - \mathbb{E}( r_b^{(a)}) \mathbb{E}( x_a ) = {\rm cov}(x_a, x_b) - \frac{{\rm cov}(x_a, x_b)}{{\rm var}(x_a)} {\rm var}(x_a) = 0$
ですから、説明変数と残差は常に独立なのではないでしょうか？

$\mathbb{E}(XY)=\mathbb{E}(X)\mathbb{E}(Y)$ は独立じゃなくて無相関の定義だよ…独立は $\mathbb{P}(X, Y)=\mathbb{P}(X)\mathbb{P}(Y)$ だから、しいてかくなら以下で、説明変数と残差は常に無相関ではあっても因果の向きを誤ると独立じゃないよ。

$\displaystyle \mathbb{P}( r_2^{(1)}, x_1 ) = \mathbb{P}( r_2^{(1)}| x_1 )\mathbb{P}(x_1) = \mathbb{P}( r_2^{(1)})p(x_1)$
$\displaystyle \mathbb{P}( r_1^{(2)}, x_2 ) = \mathbb{P}( r_1^{(2)}| x_2 )\mathbb{P}(x_2) = \mathbb{P}( r_1^{(2)}| e_2 )\mathbb{P}(e_2 | x_2)\mathbb{P}(x_2)$

だって（１）の $x_2$ を $x_1$ で説明しようとするなら説明できない分は $e_2$ で $x_1$ と独立だけど、（２）の $x_1$ を $x_2$ で説明しようとするなら説明できない分は $x_2$ と従属だよ、雑だけど。定理4.1をつかうと従属であると示せるのか…。

（１） $\displaystyle x_2 = b_{21} x_1 + e_2$
（２） $\displaystyle x_1 = \frac{x_2}{b_{21}} - \frac{e_2}{b_{21}}$

なんと…騙されました…。それで、114ページの定理4.1は、なんだかモデルが独立成分分析のようですが…どういうことですか？

定理のすぐ下にあるように対偶の方がいっていることがわかりやすいよね。「マイク1とマイク2が独立な $q$ 人の声を拾っているとする。ある $j$ さんの声が非ガウス分布にしたがい、マイク1もマイク2も $j$ さんの声を拾っている（ $\alpha_j \beta_j \neq 0$ ）ならば、マイク1とマイク2は必ず従属になる」。対偶じゃない方でいうなら、「もしマイク1とマイク2が独立なら、どちらのマイクにも声が拾われている $j$ さんの声はガウス分布にしたがっている」だね。

どちらのマイクにも声が拾われている人がいるのに独立になるんですか？？理解しがたいですね…。

証明はのってないね。まあこれを認めれば因果的順序が一番最初の変数を推定できる。独立かどうかを測るから相互情報量をつかうんだね。相互情報量の計算には結合エントロピーが要るけど、117ページのように相互情報量の差を求めにいけば結合エントロピーが消えるね。それで、因果的順序が一番最初の変数を推定できたら、その変数及びその変数の寄与を除いてまた因果的順序が一番最初の変数を探すことを繰り返すんだね。

5章以降があればつづく