2020-02-25

雑記：ダルモア・スキットビッチ定理の証明

参考文献（この記事はほとんどこの文献の内容のまま）： http://ee.sharif.edu/~bss/DarmoisTheorem.pdf

統計的因果探索手法の LiNGAM で実際に因果グラフを推測するアプローチの一つに、「正しくない因果的順序で回帰すると（原因の変数を結果の変数で回帰すると）、説明変数と残差が独立にならない ⭐」ことを利用して因果的順序を特定するものがありますが、⭐ を保証するダルモア・スキットビッチ定理の証明が気になったのでウェブ上で拾った参考文献の自分の訳をメモします。解釈の誤りは自分に帰属します。定理1と定理2の証明はありません（参考文献にも）。お気付きの点がありましたらコメントでご指摘いただけますと幸いです。

なお、ダルモア・スキットビッチ定理は以下の定理3で、「確率変数 $Y_1, Y_2$ がどちらも、非ガウス分布にしたがう $X_j$ の1次の項を含んでいる（係数 $a_j, b_j$ がどちらもゼロでない）ならば、 $Y_1$ と $Y_2$ は独立にならない」といっても同じです。LiNGAM の文脈では、 $X_j$ を結果変数の誤差項、 $Y_1$ を結果変数、 $Y_2$ を結果変数で原因変数を回帰したときの残差とすれば、 $a_j$ も $b_j$ も非ゼロであり、結果変数と誤差項が独立になりません。

補題1

$f_1, f_2, \cdots, f_N$ を何回でも微分可能な関数とし、任意の $x, y$ について以下のようにかけるとする。 $f_1(a_1x + b_1y) + f_2(a_2x + b_2y) + \cdots + f_N(a_Nx + b_Ny) = A(x) + B(y) \quad \forall x, y \tag{1}$ ただし、 $a_1, a_2, \cdots, a_N, b_1, b_2, \cdots, b_N$ はゼロでない定数であり以下を満たす。 $a_i b_j - a_j b_i \neq 0 \qquad \forall \, i \neq j$ このとき $f_1, f_2, \cdots, f_N$ は高々 $N$ 次の多項式である。

補題1の証明

$(1)$ 式において、以下のように、 $a_N x + b_N y$ は一定に保ちつつ $x$ と $y$ を変化させることを考える。
$x \leftarrow x + \delta_1^{(1)}$
$y \leftarrow y + \delta_2^{(1)}$
$a_N \delta_1^{(1)} + b_N \delta_2^{(1)} = 0$ （上の3番目の式が表す直線上の原点以外から $(\delta_1^{(1)}, \delta_2^{(1)})$ を選べばよい。）ここで、 $\epsilon_i^{(1)} = a_i \delta_1^{(1)} + b_i \delta_2^{(1)}$ は $i \neq N$ では $0$ にならない（∵ $a_i b_N - a_N b_i \neq 0$ ）。 $x$ と $y$ を変化させた後の $(1)$ 式から変化させる前の $(1)$ を引くと以下のようになる。
$\Delta_{\epsilon_1^{(1)}} f_1(a_1x + b_1y) + \Delta_{\epsilon_2^{(1)}} f_2(a_2x + b_2y) + \cdots +\Delta_{\epsilon_{N-1}^{(1)}} f_{N-1}(a_{N-1}x + b_{N-1}y) = A_1(x) + B_1(y) \quad \forall x, y$ $\Delta_h f(x)$ は $\Delta_h f(x) = f(x + h) - f(x)$ を意味する。同様の変化を $f_1$ の項だけになるまで（ $N-1$ 回目まで）繰り返すと以下のようになる。
$\Delta_{\epsilon_1^{(N-1)}} \cdots \Delta_{\epsilon_1^{(1)}} f_1(a_1x + b_1y) = A_{N-1}(x) + B_{N-1}(y) \quad \forall x, y$ $N$ 回目には $y$ を固定して $x$ だけ変化させる。 $\Delta_{\epsilon_1^{(N)}} \cdots \Delta_{\epsilon_1^{(1)}} f_1(a_1x + b_1y) = A_{N}(x) \quad \forall x, y$ $N + 1$ 回目には $x$ を固定して $y$ だけ変化させる。 $\Delta_{\epsilon_1^{(N+1)}} \cdots \Delta_{\epsilon_1^{(1)}} f_1(a_1x + b_1y) = 0 \quad \forall x, y$ ここで、各回の変化 $\epsilon_1^{(1)}, \cdots, \epsilon_1^{(N+1)}$ は好きな値にとれるが、どのような値にとっても結局 $\Delta_{\epsilon_1^{(N+1)}} \cdots \Delta_{\epsilon_1^{(1)}} f_1(a_1x + b_1y) = 0$ になる。これは $f_1$ の $N+1$ 次の導関数が $0$ ということに他ならない。したがって、 $f_1$ は高々 $N$ 次の多項式である。 $f_2, \cdots, f_N$ についても同様。

定理1（ Lévy-Cramer ）

$X_1, X_2$ を互いに独立な確率変数とし、 $Y = X_1 + X_2$ とする。このときもし $Y$ がガウス分布にしたがうならば、 $X_1$ と $X_2$ もガウス分布にしたがう。

定理1の証明は参考文献にはない。以下の2～7ページにありそう。読んでない。
[1810.01768] Three remarkable properties of the Normal distribution

定理2（ Marcinkiewics-Dugué ）

$e^{p(\omega)}$ （ $p(\omega)$ は多項式）の形の特性関数をもつ確率変数は、定数確率変数かガウス分布にしたがう確率変数だけである（したがって多項式の次数は2以下である）。

定理2の証明も参考文献にはない。同じ主張が以下の記事にもある。
Cumulants - Scholarpedia
Marcinkiewicz (1939) は探せば出てくるがフランス語なので読めない。
おそらく以下のサーベイにもこの定理が載っているがこちらはたぶんフリーアクセス版がない。
A survey of the theory of characteristic functions | Advances in Applied Probability | Cambridge Core

定理3（ Darmois-Skitovic ）

$X_1, \cdots, X_N$ を互いに独立な確率変数とし、 $Y_1, Y_2$ を以下のように定義する。 $\begin{cases} Y_1 &= a_1 X_1 + \cdots + a_N X_N\\ Y_2 &= b_1 X_1 + \cdots + b_N X_N \end{cases}$ このときもし $Y_1$ と $Y_2$ が独立ならば、 $a_i b_i \neq 0$ であるようなすべての $i$ で $X_i$ はガウス分布にしたがう。

定理3の証明

$a_i b_j - a_j b_i \neq 0 \quad \forall \, i \neq j$ としても一般性を失わない（∵ もし $a_i b_j - a_j b_i = 0$ になる $i, j$ があれば $a_i X_i + a_j X_j$ という新しい確率変数にマージすればよい；もしこのマージした確率変数がガウス分布にしたがうならば定理1よりマージする前の確率変数もガウス分布にしたがう）。
いま $(Y_1, Y_2)$ の特性関数は以下のようになる。
$\Phi_{Y_1, Y_2}(\omega_1, \omega_2) = E \bigl[e^{j(\omega_1 Y_1 + \omega_2 Y_2)} \bigr] = E \bigl[e^{j \sum_{i=1}^N (\omega_1 a_i + \omega_2 b_i) X_i} \bigr] = \prod_{i=1}^N \Phi_{X_i} (a_i \omega_1 + b_i \omega_2)$
他方、 $Y_1, Y_2$ が独立であることから、 $(Y_1, Y_2)$ の特性関数は以下のようにならなければならない。
$\Phi_{Y_1, Y_2}(\omega_1, \omega_2) = \Phi_{Y_1}(\omega_1) \Phi_{Y_2}(\omega_2)$
したがって、
$\prod_{i=1}^N \Phi_{X_i} (a_i \omega_1 + b_i \omega_2) = \Phi_{Y_1}(\omega_1) \Phi_{Y_2}(\omega_2)$
両辺の対数をとって（ $\Psi$ は第2キュムラント母関数）、
$\sum_{i=1}^N \Psi_{X_i} (a_i \omega_1 + b_i \omega_2) = \Psi_{Y_1}(\omega_1) + \Psi_{Y_2}(\omega_2)$
ここで、左辺の $a_i b_i = 0$ であるような $i$ については右辺に移項する。それで補題1を適用すると、 $a_i b_i \neq 0$ であるような $i$ で $\Psi_{X_i}$ は多項式であることがわかる。よって定理2より $X_i$ はガウス分布にしたがう（分散がゼロでないと仮定しないならば定数確率変数でもありうると思う）。

2020-02-21

統計的因果探索：ノート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章以降があればつづく

2020-02-20

統計的因果探索：ノート1章、2章、3章

本読み統計的因果探索

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

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

著者のスライド

まとめ

1章：【悲報】因果効果、相関係数からわからない
2章：構造的因果モデルの導入、それを用いた平均因果効果の定義、平均因果効果が推定できるランダム化実験
3章：ランダム化実験できないときはどうやって平均因果効果を推定するの？ → 因果グラフさえ特定すれば勝てる（そのノードXのノードYに対する平均因果効果はノードXとノードXの親たちでノードYを回帰した偏回帰係数だから）→ じゃあどうやって因果グラフ推定するの？ → 構造的因果モデルに線形関数・非ガウスノイズを仮定すれば観測変数の分布から因果グラフを識別できて勝てる（LiNGAM）（非線形関数を許す手法もある；p56）（逆に関数とノイズに何も仮定しなかったり、線形関数・ガウスノイズだと因果グラフを識別できない）

Chapter1 の例（13ページの図1.4）は、同じ相関係数なのに色々な因果のパターンがあるという例だけど、 $\rho=0.79$ でありうるすべてのパターンを図示すると以下のグラフの青線になるかな（もちろん Y → X の向きも同様ね）。

未観測共通原因 $Z$ がないなら素直に $X$ から $Y$ に 0.79 の因果効果があるけど、 $X$ も $Y$ も未観測共通原因 $Z$ から $a$ の因果効果を受けているときは、そのことが相関係数に $a^2$ を上乗せしてしまう。だから、 $X$ の $Y$ に対する因果効果は見かけの相関係数からそれを差し引かないといけない（ $b = \rho - a^2$ ）。だから、 $b$ は 0.79 より小さいかもしれないし、0（因果関係なし）かもしれないし、下手したら負かもしれない。

見かけ上は正の相関があるが本当の因果効果は負である、なんてことがあるでしょうか。

そうだな、現実とは違うかもしれないけど、「ゲームソフトをたくさんもっている児童ほどテストの成績がよい」という正の相関があったとする。でも、実は「ゲームソフトをたくさんもっている」にも「テストの成績がよい」にも「家が裕福である」という交絡因子があまりにも大きく効いていたとする。そしたら、本当は「ゲームソフトをどれだけもっているか」の「テストの成績」に対する効果が負であったとしても正の相関が観測されるかもしれない。あとは、前の記事にもかいたけど、「若手の医師ほど手術の成功率が高い」とか。ベテランの医師ほど難しい手術をまかされるならね。…あ、でも、後者の場合は下図の左ではなくて右なのかな。「X: 経験年数」「Y: 手術の成功率」「Z: まかされる手術の難易度」だとしたら、X → Z だと思うし。

f:id:cookie-box:20200220110030p:plain:w360

なるほど…Chapter1 を要約すると「因果効果の大きさ、わからない」といったところでしょうか。

それで Chapter2 に入ると、構造的因果モデルの流儀で「この集団において、 $x$ は $y$ の原因になる」とは、 $p(y|{\rm do}(x=d)) \neq p(y|{\rm do}(x=c))$ を満たす $c$ と $d$ が存在する。ということで、平均因果効果は $\mathbb{E}(y|{\rm do}(x=d)) - \mathbb{E}(y|{\rm do}(x=c))$ ですね。他方、「個体Aにおいて、 $x$ は $y$ の原因になる」は以下ということです。 $y^{\rm (A)}_{x=d} \neq y^{\rm (A)}_{x=c}$ を満たす $c$ と $d$ が存在する。

それで、(2.3) (2.4) 式の構造方程式モデルで $x$ の値を $0$ から $1$ に変化させたときの平均因果効果は $\mathbb{E}(y|{\rm do}(x=1)) - \mathbb{E}(y|{\rm do}(x=0)) = \mathbb{E}(f_y(1, e_y)) - \mathbb{E}(f_y(0, e_y))$ ですが、 $f_y(0, e_y)$ の分布と $f_y(1, e_y)$ の分布を同時に観測することはできないのですよね（22ページ）。そこで、 $x$ をランダムに $0$ か $1$ に決めることを考えます（ベルヌーイ分布で；テキストでは成功確率 1/2 とされていますが、ここの文脈に限っては 1/2 である必要もないはずです）。このとき、 $f_y(0, e_y), \, f_y(1, e_y)$ のしたがう分布を $x$ で条件付けても分布の形は変わりませんね。そうすると以下のトリックがつかえます。 $\begin{split}\mathbb{E}(y|{\rm do}(x=1)) - \mathbb{E}(y|{\rm do}(x=0)) &= \mathbb{E}(f_y(1, e_y)) - \mathbb{E}(f_y(0, e_y)) \\ &= \mathbb{E}(f_y(1, e_y)|x=1) - \mathbb{E}(f_y(0, e_y)|x=0) \\ &= \mathbb{E}(f_y(x, e_y)|x=1) - \mathbb{E}(f_y(x, e_y)|x=0) \\ &= \mathbb{E}(y|x=1) - \mathbb{E}(y|x=0) \end{split}$ この式変形を日本語で説明するなら、さしずめ以下でしょう。

ランダム化実験のストーリー

ある薬がある病気の人たちを治す効果を測りたいのですが、「全員に薬を飲ませたときの病気の人の割合」と「全員に薬を飲ませなかったときの病気の人の割合」を両方観測することはできません。前者を観測した後時間を巻き戻して後者を観測できればいいのですが、時間は巻き戻せませんので。
全員にあらかじめコインを振ってもらうことにしましょう。まあ、あらかじめコインを振ってもらったところで、「全員に薬を飲ませたときの病気の人の割合」と「全員に薬を飲ませなかったときの病気の人の割合」を両方観測することはできません。
ところでこのとき「全員に薬を飲ませたときの病気の人の割合」と「全員に薬を飲ませたときの、コインが表の人の中での病気の人の割合」は等しくなるはずですね。コインを振った結果が病気が治るかどうかに影響を与えることはありませんから。同様に、「全員に薬を飲ませなかったときの病気の人の割合」と「全員に薬を飲ませなかったときの、コインが裏の人の中での病気の人の割合」も等しくなるはずです（上式の2行目）。
さらに、「全員に薬を飲ませたときの、コインが表の人の中での病気の人の割合」は、「コインが表の人に薬を飲ませてコインが裏の人に薬を飲ませなかったときの、コインが表の人の中での病気の人の割合」と言い換えられます。「全員に薬を飲ませなかったときの、コインが裏の人の中での病気の人の割合」は「コインが表の人に薬を飲ませてコインが裏の人に薬を飲ませなかったときの、コインが裏の人の中での病気の人の割合」と言い換えられます（上式の3行目）
「～ときの」がそろったのに気付くでしょう。もはや時間を巻き戻す必要がなくなったのです。

ここで仮定として、「病気が治っているかどうか」から「薬を飲むかどうか」の向きに因果があってはならないです。ここでは、後者は前者に時間的に先行しているので原因たりえないとありますが。

ときに副部長、上のランダム化実験の説明をかいていたら因果推論へのオールマイティなアプローチを思い付きました。タイムマシンの開発です！

…部長だけがこっそり過去の行動を変えることができるタイムマシンがあるんだったらもう因果推論なんかやらずに20年前に戻ってヤフー株買えばいいじゃないか。

それもそうですね。Chapter2 は構造的因果モデルの導入でした。Chapter3 は、この本のタイトルである「統計的因果探索」の基礎ということですが、そもそも統計的因果探索とは因果グラフを未知としてどうすれば因果グラフが推測できるかを取り扱う研究なのですね。確かに、「睡眠障害」と「抑うつ気分」はどちらが原因でどちらが結果なのか、はたまた因果関係はないのかわかりませんね。そして49ページにきて、Chapter1 の例（3つの因果グラフ）に戻ってきました。データが3つの因果グラフのうちどれから生成されたかを推測することを「因果探索の基本問題」とよぶと。そして「因果探索の基本問題」はノンパラメトリックアプローチでは解けない？では、パラメトリックに、かつ関数には線形性を・外生変数の分布にはガウス分布を仮定したら、やっぱり解けないと…。

13ページの図1.4の3つの観測変数の分布は全部同じで区別できないってことだね。

それで、このテキストで主に取り扱うのはセミパラメトリックに関数形のみを仮定し外生変数の分布には仮定をおかないアプローチで、特に関数に線形性・外生変数の分布に非ガウスを仮定する LiNGAM (Linear Non-Gaussian Acyclic Model) を解説すると。これで「因果探索の基本問題」が解けるということですが。「ガウス分布でさえなければ、どんな連続分布でもかまいません（53ページ）」というのもすごいですね。逆にいうと、それほどまでにガウス分布はデータがどうやって生まれたかを包み隠してしまうのでしょうか…。

それで、58ページからは未観測共通原因 $Z$ が観測済みであると仮定していますね。そうすると $X, Y, Z$ の因果グラフの候補は27通り（双方向の因果関係は考えないので；46ページ）…多いですね…非巡回の仮定を追加しても25通り…全然減ってないじゃないですか…。

しかし、未観測共通原因がなく因果関係が非巡回であれば、 $x_i = f_i ({\rm pa}(x_i), e_i)$ (3.2) と、自分の親と自分の誤差変数から影響を受けるのみですね。さらに線形関数を仮定すれば $x_i = \sum_{x_j \in {\rm pa}(x_i)} b_{ij} x_j + e_i$ (3.3) です。これをさらに行列で $x = Bx + e$ (3.5) とかけば行列 $B \in \mathbb{R}^{p \times p}$ の1つのゼロ・非ゼロのパターンが1つの有効非巡回グラフに対応すると。それはそうですね。そして、ここでの統計的因果探索の目的が以下のように定式化されます。データ行列 $X \in \mathbb{R}^{p \times n}$ から構造方程式の係数行列 $B \in \mathbb{R}^{p \times p}$ を推定する。この $B$ が推定できればどの変数とどの変数が親子関係にあるかわかります。そうなれば、 $x_i$ に介入したときの $x_j$ の期待値が以下で求めることができるとありますね（ただし $x_j$ は $x_i$ の祖先ではないとします）。 $\mathbb{E}(x_j | {\rm do}(x_i = c)) = \mathbb{E}_{{\rm pa}(x_i)} [ \mathbb{E}(x_j | x_i = c, \, {\rm pa}(x_i)) ]$ えっと、右辺、なぜこうなっているのでしょう？一度 $x_i = c$ と $x_i$ の親で条件付けて $x_j$ の期待値をとってから、 $x_i$ の親について期待値をとるんですか？

単純に $\mathbb{E}(x_j | x_i = c)$ ではいけないのはわかるよね。これでは「 $x_i$ の親たちが $x_i = c$ にしてくれるようなとき」について期待値をとっていることになる。でもこれだと、 $x_i$ の親の誰かから $x_j$ に対する因果関係はあるが実は $x_i$ と $x_j$ の間に因果関係はないというときにも、 $x_i = c$ の $c$ に応じて $x_j$ の期待値は変わってしまう。でも、こんなのは $x_i$ の $x_j$ に対する因果効果じゃない。 $x_i$ を変化させたって $x_j$ は変化しないからね。本当に $x_i$ を変化させたときの $x_j$ の変化を知りたいなら、「 $x_i$ さんの親たちの実現値にかかわらず $x_i$ を $c$ にしたときの」という条件付けが必要だ。まず $x_i$ さんの親たちのある実現値に対して $x_j$ の期待値を計算したのが $\mathbb{E}(x_j | x_i = c, \, {\rm pa}(x_i))$ だよね。他方、 $x_i$ さんの親は $x_i$ を経由せずに $x_j$ に因果効果をもたらしている可能性はある。そこは期待値計算のときに考慮しなきゃだから、そこは普通に期待値をとって結局 $\mathbb{E}_{{\rm pa}(x_i)} [ \mathbb{E}(x_j | x_i = c, \, {\rm pa}(x_i)) ]$ になる。

む、そういわれるとそうですね…そして、 $x_i$ から $x_j$ への平均因果効果は $\mathbb{E}(x_j | {\rm do}(x_i = d))$ と $\mathbb{E}(x_j | {\rm do}(x_i = c))$ の差であり、結局 $\alpha_{ji}(d - c)$ になるんですね。この $\alpha_{ji}$ は $x_j$ を $x_i,\, {\rm pa}(x_i)$ で回帰したときの $x_i$ の偏回帰係数です。

まとめると、「 $x_i$ から $x_j$ への平均因果効果」を求める手順としては以下になるのかな。

構造方程式の係数行列 $B$ を推測し、因果グラフを描く。
$x_j$ が $x_i$ の祖先ではないことを確認する（祖先であったら因果の向きに反する）。
$x_i$ の親たち ${\rm pa}(x_i)$ をリストアップする。
$x_j$ を $x_i,\, {\rm pa}(x_i)$ で回帰する。
$x_i$ の偏回帰係数 $\alpha_{ji}$ が $x_i$ から $x_j$ への平均因果効果である。

72ページの図3.11には逆に説明変数に加えてはいけない例というのがありますが、その手順にきちんとしたがうなら誤って含めてしまうことはないですよね…。

74ページからはノンパラメトリック、パラメトリック、セミパラメトリックアプローチの比較ですね。というかまず、因果グラフが有効非巡回グラフならば因果的マルコフ条件が成り立つと。そして、忠実性？「変数間の条件付き独立性が因果的マルコフ条件から導かれるものだけである」？それはそうなんじゃないですか？因果的マルコフ条件から導かれない条件付き独立性って何ですか？？

77ページの図3.13、因果グラフだけみれば因果的マルコフ条件から導ける条件付き独立性がないのに、実は $p(x,y,z)=p(z|y)p(y|x)p(x)$ になってるね…これは… $x$ から送り出された $x$ と $-x$ が $z$ で合流して起きた負の奇跡かな？実質的に $x$ から $z$ への有向辺がないのと同じになってる。

なるほど？しかし、それを仮定する必要があるんですか？だって、その77ページの図って、まさに副部長がいった通り「実質的に $x$ から $z$ への有向辺がないのと同じ」なんでしょう？じゃあ $x$ から $z$ への有向辺がない因果グラフを推測したって問題ないじゃないですか。それとも、誤差の議論などで不具合が出てくるんでしょうか。確かにあるはずの枝を見落とすと平均因果効果がくるってきそうですが…まあいいです。それで、各アプローチで因果グラフを特定できるかどうかを概観すると以下のようになるんですね。

ノンパラメトリックアプローチ： マルコフ同値類までしか特定できないので因果グラフを一意に推測できないことがある。以下の手順による。

観測の変数の分布から条件付き独立性を見出す。
その条件付き独立性を与える因果グラフを探索する。

PCアルゴリズム
GESアルゴリズム

（線形性・ガウス性を仮定する）パラメトリックアプローチ： 同じ観測変数の分布を与える因果グラフを識別できない。

（平均因果効果のとりうる範囲を推定しようとする試みなどはある。）

セミパラメトリックアプローチ（LiNGAM）： 誤差変数が非ガウスなら因果グラフを一意に特定できる。詳細は次章以降。

各因果グラフの場合の観測変数の分布を求める（どうやって誤差分布を特定する？？）。
実際の観測変数の分布と比較する。

上でパラメトリックアプローチは線形・ガウシアンという制約を課されていますが、まあガウシアンでない場合はすべてセミパラメトリックアプローチに吸収されるのだとして、関数が非線形、ノイズがガウシアンなどであったらどうなるんでしょうか…まあ面倒そうなのでいいですが…。

4章以降があればつづく

2020-02-16

論文読みメモ： Temporal regularized matrix factorization for high-dimensional time series prediction（その1）

論文読み

以下の論文を読みます。

Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon. Temporal regularized matrix factorization for high-dimensional time series prediction. In Advances in neural information processing systems, pages 847–855, 2016. https://papers.nips.cc/paper/6160-temporal-regularized-matrix-factorization-for-high-dimensional-time-series-prediction

※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。

前回：論文読みメモ： Think Globally, Act Locally: A Deep Neural Network Approach to High-Dimensional Time Series Forecasting（その1） - クッキーの日記

前回の論文の引用文献の話が長くなりそうだったから記事をスピンアウトするね。つまり、高次元時系列予測モデルの先行研究であって、ニューラルネットを利用するものであって、グローバルなパターンも学習するものの2つ目の方ね。論文の図のまんまだけど、時系列を下図のように $Y \approx FX$ と行列分解（MF）するんだね。 $F$ の各行は $i$ 番目の商品の特徴で、 $X$ の各列はその時刻の特徴って感じだね。

f:id:cookie-box:20200216104950p:plain:w480

こんな行列分解を求めるには以下を解くんだけど、2番目の項と3番目の項は正則化項だね。 $\displaystyle \underset{F, X}{\rm min} \sum_{(i, t) \in \Omega} (Y_{it} - f_i^\top x_t)^2 + \lambda_f \mathcal{R}_f(F) + \lambda_x \mathcal{R}_x(X)$ でも通常（とは）の行列分解みたいに $\mathcal{R}_x(X) = \| X \| _F$ （ $\| \cdot \| _F$ は後で出てくるけどフロベニウスノルムだね…最初 $Y \approx FX$ の $F$ と関係ある何らかのノルムなのかと思った…）なる正則化は適切ではないといっているね。もし $Y$ がアイテム×ユーザという行列だったら、各ユーザがどのアイテムをより好きかを表現できればいいけど、いまはアイテム×時刻という行列だから、「各時刻にどのアイテムがより強い（？）」だけでは不十分で、次の時刻にその強さがどう変わるかも必要だもんね。

多変量時系列分析に行列分解を用いた先行研究ではこういう時間方向の依存性をグラフベースで取り扱ってラプラス正則化を適用しているみたいだけど、これだと2時点間の負の相関を考慮できなかったり、そもそも時間依存性が明示的に手に入っていない場合は推測しないといけなかったり、欠損値には強いものの予測は不得手だったりするって。

でも、この論文で提案する temporal regularized matrix factorization framework (TRMF) は「時間的正則化項」を導入することでデータドリブンで時間依存性を学習できて予測性能もあると。MF の特徴であるスケーラビリティと欠損値への強さはそのままに。あと先行のグラフベースアプローチとの統一的な見方も与えると。そうすることで既存のソルバーが活用できるって。

MFを用いる先行手法にはかなり制約があるのですね…？本論文の2節には先行手法とその限界についてまとめてありますね。

まず、AR や DLM は計算量が時系列の本数や参照ステップ数にスケールしないと。
- L次のARモデルって $O(Ln^2)$ のパラメータの推測に $O(TL^2n^4 + L^3n^6)$ かかるんですか？
- DLMってパラメータを更新するのに $O(kn^2T + k^3 T)$ かかるんですか（ $k$ は隠れ変数の次元数でしばしば $n$ よりも大きい）？なので R の DLM で数十を超える次元数を取り扱うことはできないと。
次に、MF を用いる手法ですね。まず、グラフベースの正則化をする手法が以下ですよね。3番目の文献は5番目の文献の extended version にみえるのですがどう違うのですかね…。
ただ以下の先行手法は例外で、正則化項をカスタマイズしているということですが、これは音声認識の論文でしょうか。
- Regularized Non-negative Matrix Factorization with Temporal Dependencies for Speech Denoising (2008)
まあそれで、グラフベースの正則化とは以下のようなものであるということですね。
$\displaystyle \mathcal{R}_x(X) = \mathcal{G}(X|G, \eta) := \frac{1}{2}\sum_{t \sim s}G_{ts} \| x_t - x_s \|^2 + \frac{\eta}{2}\sum_t \|x_t\|^2$ が時刻と時刻の間の枝の重みということでしょうか。は枝があるところについての和ですね。なるほど、つまり、「結びつきが強い枝でつながっているとほど同じ値であれ」という最適化ですね。2項目はよくあるように、最適な点があることを保証するためのものっぽいですね。そして、でなければなりません。それはそうでなければ最小化が実行できませんものね…。そして、実際には相関があると信ずるタイムラグの集合、例えばなどに対して適当な時不変な重みを設定して上の最適化をするのですね。しかしそれでは、負の相関を考慮できないし、そもそも明示的に時間依存構造が必要だし、予測性能もよくないのですね。
じゃあ、時間依存構造 $w_l$ そのものもデータから学びたいということになりますよね。素朴には、以下のような最適化をしたいということになります。
$\displaystyle \underset{F, X, w \geqq 0}{\rm min} \sum_{(i, t) \in \Omega} (Y_{it} - f_i^\top x_t)^2 + \lambda_f \mathcal{R}_f(F) + \frac{\lambda_x}{2} \sum_{l \in \mathcal{L}} \sum_{t:t-l>0} w_l (x_t - x_{t-l})^2 + \frac{\lambda_x \eta}{2} \sum_t \|x_t\|^2$ しかし、これを解くと任意の $l$ に $w_l=0$ となるのは明白だと…当然ですね…。なのでこのアプローチは取れません…いったいどうすれば…。

もう少しつづく

2020-02-07

論文読みメモ： Think Globally, Act Locally: A Deep Neural Network Approach to High-Dimensional Time Series Forecasting（その1）

論文読み

以下の論文を読みます。

Rajat Sen, Hsiang-Fu Yu, Inderjit S. Dhillon. Think Globally, Act Locally: A Deep Neural Network Approach to High-Dimensional Time Series Forecasting. In Advances in Neural Information Processing Systems 32, 2019.
https://papers.nips.cc/paper/8730-think-globally-act-locally-a-deep-neural-network-approach-to-high-dimensional-time-series-forecasting

※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらご指摘ください。

次回：まだ

近年の時系列データには、一系列が一個人に対応し、何千もの系列が互いに相関するようなものがありうるとありますね。なのでそれらの系列からグローバルなパターンを搾り取った上で、個別データに対するキャリブレーションをするとよい予測ができるのではないかといっています。そこでこの論文で提案するモデル、DeepGLO は「グローバルな」Matrix Factorizartion モデルと、個別の系列の性質及びそれらの相関を捉えるネットワークを組み合わせるということです。

冒頭で時系列予測は様々な分野で重要な課題だといって挙げられている例は、1番上のは本論文と同じ Amazon の方々の論文で、2番目は割と古いですね。3番目も結構古い本です。

Bayesian Intermittent Demand Forecasting for Large Inventories（retail demand forecasting）
Financial time series forecasting using support vector machines - ScienceDirect（financial predictions）
Amazon | Time-Series Forecasting | Chris Chatfield | Applied（ predicting traffic or weather patterns）

最近のデータセットは何千もの相互に相関する時系列 $Y \in \mathbb{R}^{n \times t}$ が含まれているということですが、ECサイトでいう「各カテゴリの需要」や「各商品の需要」というのを指しているのでしょうか。需要が高まりそうな商品については契約している業者に在庫の補充を促したりするのですかね？

それで、時系列予測は伝統的には個々の系列や少数の系列たちに対して AR、ARIMA、指数平滑法、Box-Jenkins法、状態空間モデルを適用する方法がとられてきたとありますね。このうち、状態空間モデルについては複数の系列を一緒に取り扱えますよね。しかし、膨大な数の系列にはスケールしないということですか。まあそんな気はしますが…。そして、すべての系列が共通にもっているパターンを上手く利用できないともあります。これは状態空間モデルでもそうなんでしょうか。よくわかりません。それらに対してニューラルネット手法として RNN、LSTM、WaveNet の話があって、これらは大規模な時系列も学習できるものの、2つの欠点があると指摘していますね。

個々の時系列の規模が異なると上手くいかない。例えば個々の製品の需要予測の場合、人気商品の需要の規模はニッチな商品のそれよりも大きい。このような場合、上記のようなニューラルネット手法では正規化が必要になる。しかし、正規化のパラメータ調整は難しい。
モデルは予測時に個別のデータの過去データしか利用していない。しかし、データはすべての時系列にわたるグローバルなパターンをもっている。

なので、これに対処したモデルも提案されている（後者は本論文と同じ方の論文ですが、このときは所属が Amazon ではないのですね）（前者はコードがありますが Python2.7 だと…）。

[1703.07015] Modeling Long- and Short-Term Temporal Patterns with Deep Neural Networks GitHub - laiguokun/LSTNet	Long- and Short-term Time-series network （LSTNet）； 2次元畳込みと再帰結合を組み合わせて、予測時に複数時系列のグローバルな情報を利用できるようにした。しかし、何千次元にはスケールしない。
https://papers.nips.cc/paper/6160-temporal-regularized-matrix-factorization-for-high-dimensional-time-series-prediction	Temporal Regularized Matrix Factorization （TRMF）；すべての時系列を「基底時系列」たちの線形結合で表現する。この「基底時系列」がグローバルなパターンを捉える。しかし、線形な時間依存性しか表現できない（？）。また、行列分解の際に近似誤差が生じ、個々のパターンの情報が失われやすい。

…グローバルなパターンを活用しようとしたモデルについてこう簡単にかかれていますが、具体的にどのようにグローバルなパターンを捉えているのですかね？

1つ目の方は論文内の図がわかりやすいけど以下のような感じだね。まず2次元畳込みをして、その出力を普通にGRUする層と、スキップしてGRUする層もつくって（例えば1時間毎といった周期性があるデータなら1時間スキップして1時間毎のレベルの変化を学習するみたいな）、それらの出力を全結合してニューラルネットの予測結果にする。それを、通常のARモデルの結果と足し合わせるみたい。

f:id:cookie-box:20200208002515p:plain:w480

2つ目はスケーラブルで欠損値があっても大丈夫らしい。Wal-mart E-commerce datasets でよい予測ができたと。

もう少しつづく