因果指標と偏正準相関分析
8
0
0
全文
(2) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. 3.1 Transfer Entropy と Granger Causality. た下での条件付きエントロピーとして表現できる．. ∫∫∫. (l). ここでは TE と GC の関係性を述べる [4]．X, Y が連. (k). p(yt , yt−1 , xt−1 ) ×. Tx→y =. 続値データであるとする．時系列データ Y が式 (7) の. (k) (l) p(xt−1 , yt |yt−1 ) (l) (k) log2 dyt dyt−1 dxt−1 (k) (l) (l) p(xt−1 |yt−1 )p(yt |yt−1 ). = Hx(k) |y(l) + Hyt |y(l) − Hx(k) ⊗yt |y(l) .(5) t−1. t−1. t−1. t−1. t−1. 2.2 Continuous Transfer Entropy. AR モデルに従うと仮定した場合，yt の条件付き確率分布 (l). p(yt | yt−1 ) は以下の式で与えられる． ( )2 T (l) y − a y t t−1 1 √ exp − . 2 2 2σ 2πσy y. (9). 時系列データがガウス分布に従うと仮定した場合の TE (l). （Tx→y ）は次のように表される．. ¯ ¯¯ ¯ ¯Σ ¯¯ ¯ (l) (k) (l) (k) ¯ ¯Σ (l) (l) ¯ {yt−1 ,xt−1 }{yt−1 ,xt−1 } {yt ,yt−1 }{yt ,yt−1 } ¯ 1 ¯¯ ¯ , log2 ¯ ¯Σ ¯¯ ¯ 2 (l) (k) (l) (k) ¯ ¯Σ (l) (l) ¯ ¯ {yt ,yt−1 ,xt−1 }{yt ,yt−1 ,xt−1 } yt−1 yt−1 (6). ( ここで，Σy(l). (l). (l). ). は確率分布 p yt−1 の共分散行列で，. ¯ ¯ ¯ ¯ ¯Σ¯ はその行列式を意味する．これは Continuous Transfer Entropy（以下 CTE と略す）と呼ばれている [2]． t−1 yt−1. 2.3 Granger Causality Granger Causality（以下 GC と略す）[3] は，時系列デー タの予測を行う際の，自身のデータのみを用いた自己回帰 モデルと，自身のデータの他にもう 1 つ別のデータを加え た混合回帰モデルとの精度を比較する手法である．. yt = a0 +. +. (y) ϵt ,. (l). 布となる．いま，yt−1 の平均を µy(l) ，分散共分散行列を t−1. (l). Σy(l). (l) t−1 yt−1. とすると，同時確率分布 p(yt , yt−1 ) は以下のよ. うな多次元正規分布となる．. 1 √ ¯ ¯ l+1 (2π) σy2 ¯Σy(l). (l) t−1 yt−1. ( ) 1 (l+1) T (l+1) Xψ exp − ψ t ¯ 2 t ¯ ¯ (10). T y − a µ (l) t yt−1 (l+1) , = (l) where ψt yt−1 − µy(l) t−1 1 aT − 2 2 σ σy . X = ya aaT −1 − σ2 + Σ 2 (l) (l) σ y. (11). (12). yt−1 yt−1. y. ¯ ¯. この分布の分散共分散行列の行列式は σy2 ¯Σy(l). まず，自己回帰モデル（AR モデル）を考える． (l) aT yt−1. また，AR モデルでは yt−1 の確率分布は多次元正規分. (l) t−1 yt−1. ¯ ¯ ¯ とな. る．よって，以下の式が成り立つ．. (7). ここで，a0 は定数項，a ∈ R は最小二乗法により算出さ. ¯ ¯ ¯ ¯ ¯Σ ¯ ¯ 2¯ Σ (l) (l) (l) (l) ¯ . = σ ¯ {yt ,yt−1 ¯ ¯ y }{yt ,yt−1 } yt−1 yt−1. (13). l. (y). れる回帰係数である．また，ϵt. は平均 0，分散 σy2 の正規. 分布に従う誤差である． それとは別に次のような混合回帰モデルを考える． (l). (k). (y|x). yt = b0 + bT yt−1 + cT xt−1 + ϵt. .. (8). AR モデルと同様，b0 は定数項, b ∈ Rl , c ∈ Rk は最小二 (y|x). 乗法により算出される回帰係数である．また，ϵt. は平. 2 均 0，分散 σy|x の正規分布に従う誤差である．. (l). (k). {yt , yt−1 , xt−1 } の同時確率分布についても，式 (8) の混 合回帰モデルを仮定すると，以下の式が成り立つ．. ¯ ¯ ¯Σ ¯ (l) (k) (l) (k) ¯ = ¯ {yt ,yt−1 ,xt−1 }{yt ,yt−1 ,xt−1 } ¯ ¯ ¯ 2 ¯ (l) (k) (l) (k) σy|x ¯Σ{yt−1 ,xt−1 }{yt−1 ,xt−1 } ¯ .. 回帰モデルに従う時系列データの同時確率分布は多次元 正規分布になることから，式 (13) および式 (14) を CTE の 定義式である式 (6) に代入すると，以下の式が導出される．. 回帰モデルのモデリング精度は, 誤差の分散によって評 2 価できる．つまり，σy2 と σy|x を比較することで X から Y. への因果性を計ることができる．Y のモデリング精度が. X を回帰モデルに取り込むことによって向上し，分散が小 (σy2. Tx→y =. σy2 1 log2 2 . 2 σy|x. 算出した指標となっている．. 用いられている．. 3.2 Transfer Entropy と偏正準相関分析. 本章では因果指標の関係性を統一的な解釈を与，因果指 標の拡張を考える．. c 2014 Information Processing Society of Japan ⃝. (15). この指標は 2 つの回帰モデルのモデリング精度を比較する. 2 さくなった > σy|x ) ならば，X は Y に影響を持つと 2 2 2 言える．しかし，σy と σy|x との比較方法は σy2 − σy|x や 2 2 σy|x /σy など複数の方法が考えられ，実際に様々な方法が. 3. 因果指標の関係性と拡張. (14). という枠組みの中で，X が Y に与えている平均情報量を. ここでは，TE が偏正準相関分析に帰着されることを 示す [5]．確率変数 x, y, z に対して，偏相関係数の考え方 と同様にして z の影響を除いた x と y で正準相関分析を 行うのが偏正準相関分析（Partial Canonical Correlation. 2.
(3) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. Analysis，以下 PCCA と略す）[6] である．PCCA は以下. ここで λ はラグランジュの未定乗数である．a と λ に関し. の一般化固有値問題を解くことで求められる. て，式 (23) の停留条件により以下の式が導かれる．. Σxy|z Σ−1 yy|z Σyx|z a = λΣxx|z a, Σyx|z Σ−1 xx|z Σxy|z b. Σyt x(k) |y(l) Σ−1 (k). (16). t−1. = λΣyy|z b.. t−1. (k). (l). xt−1 xt−1 |yt−1. Σx(k) yt |y(l) a t−1. (17). t−1. = λΣyt yt |y(l) a. (24) t−1. ここで，Σuv|w = Σuv −. Σuw Σ−1 ww Σwv. であり，Σuv は u. これは一般化固有値問題であり，PCCA と等価となる．こ. と v との共分散行列である．式 (16) と (17) は共通の固有. の問題を解くことによって，Dy 個の固有値と固有値に対. 値を持つ．. y ，(λ1 ≥ · · · ≥ λDy ≥ 0) を 応する固有ベクトル {λi , ai }i=1. D. まず，2 つの時系列の確率変数 x と y の次元をそれぞれ. 得る．よって TE は次のように表される．. Dx と Dy する．また 2 つの変数ともにガウス分布に従う とする．この時，各ガウス分布のエントロピーは以下のよ うに表現される．. Tx→y =. kDx log2 (2πe) 2 ¯ ¯ 1 ¯ ¯ + log2 ¯Σx(k) x(k) |y(l) ¯ , (18) t−1 t−1 t−1 2 ¯ ¯ Dy 1 ¯ ¯ Hyt |y(l) = log2 (2πe) + log2 ¯Σyt yt |y(l) ¯ ,(19) t−1 t−1 2 2 kDx + Dy Hx(k) ⊗yt |y(l) = log2 (2πe) t−1 t−1 2 ¯ ¯ ¯ ¯ 1 ¯ ¯¯ ¯ + log2 ¯Σx(k) x(k) |y(l) ¯ ¯Σyt yt |{x(k) ,y(l) } ¯ ,(20) t−1 t−1 t−1 t−1 t−1 2 こ こ で ，Σuu|{v,w} = Σuu|w −. Σuu|v − Σuw|v Σ−1 ww|v Σwu|v. (25). 値を評価することと等価となる． (k). (k). また，xt−1 の線形写像 bT xt−1 も考えられる．この場 合，bT Σx(k) x(k) |y(l) b = 1 の拘束条件を用いることで次の t−1. t−1. t−1. 式が導かれる．. =. Σx(k) yt |y(l) Σ−1 t−1. t−1. (l). yt yt |yt−1. Σyt x(k) |y(l) b t−1. t−1. = λΣx(k) x(k) |y(l) b.(26) t−1. t−1. t−1. 式 (24) と (26) は式 (16) と (17) と同様に，共通の固有値 を持つ．そのために，行列の次元や要求される固有ベクト. 1 ¯ = log2 ¯ ¯Σ ¯ 2 (l) ¯ yt yt |{x(k) ¯ t−1 ,yt−1 } ¯ ¯ ¯Σ ¯ (k) (k) (l) ¯ ¯ xt−1 xt−1 |yt−1 1 ¯. = log2 ¯ ¯Σ ¯ 2 ¯ x(k) x(k) |{yt ,y(l) } ¯ t−1. 2. 1 , 1 − λi. よって，ガウス分布の仮定の下では，TE は PCCA の固有. である．式 (18), (19), と (20). ¯ ¯ ¯Σ ¯ (l) ¯ ¯ yt yt |yt−1. t−1. i=1. log2. ここで Tx→y は i 番目の部分空間における TE の値である．. を式 (5) に代入すると，以下の式を得る．. Tx→y. Dy ∑ 1. {i}. t−1. Σuv|w Σ−1 vv|w Σvu|w. {i} = Tx→y. i=1. Hx(k) |y(l) = t−1. Dy ∑. (21). ルによってどちらの問題を解くのかを選択することができ る．固有ベクトル a と b は次の関係で表現される．. t−1. T. ここで，st = a yt つまり yt の線形写像に関して TE を 最大化することを考える．この射影は多次元データから情. 1 ai = √ Σ−1 (l) Σyt x(k) |y(l) bi , t−1 t−1 λ yt yt |yt−1 1 −1 bi = √ Σ (k) (k) (l) Σx(k) yt |y(l) ai . t−1 λ xt−1 xt−1 |yt−1 t−1 3.3 Transfer Entropy と相互情報量. 報の流れを多く含む線形部分空間を求めることに対応す. TE と時間遅れのある相互情報量（Time Delayed Mutual. る．式 (21) を用いることで，x から st への TE は次のよ. Information 以下 TDMI と略す）をシステマティックに評. うに書くことができる．. 価する [5]．このために，TDMI を xt−k の代わりに xt−1. Tx→st. (k). aT Σyt yt |y(l) a 1 1 1 t−1 = log2 T = log2 , 2 a Σyt yt |{x(k) ,y(l) } a 2 1−ρ t−1. t−1. を用いて拡張する．つまり，. ∫∫. Mx(k) yt = t−1. ρ=. aT Σyt x(k) |y(l) Σ−1 (k) t−1. t−1. (k). (l). xt−1 xt−1 |yt−1. Σx(k) yt |y(l) a t−1. aT Σyt yt |y(l) a. t−1. (k). p(xt−1 , yt ) × (k). log2. . (22). t−1. Tx→st と ρ の最大化は同じであることに注意されたい．こ こで aT Σy. (l) t yt |yt−1. a = 1 の拘束条件のもとで式 (22) の分. 子の最大化を考える．これに対応するラグランジュ関数は 以下のように与えられる．. J = aT Σyt x(k) |y(l) Σ−1 (l) a (k) (k) (l) Σ (k) xt−1 xt−1 |yt−1 xt−1 yt |yt−1 t−1 t−1 ( ) +λ 1 − aT Σyt yt |y(l) a , (23) t−1. c 2014 Information Processing Society of Japan ⃝. 1 log2 2. p(xt−1 , yt ) (k) p(xt−1 )p(yt ). ¯ ¯ ¯ ¯ ¯Σ yt yt ¯. ¯ ¯ ¯Σ ¯ ¯ yt yt |x(k) ¯ t−1 ¯ ¯ ¯Σ ¯ (k) (k) ¯ ¯ xt−1 xt−1 1 ¯ ¯. = log2 ¯Σ ¯ 2 (k) ¯ x(k) ¯ t−1 xt−1 |yt =. (k). dxt−1 dyt. (27). TE と PCCA との関係と同様に， st = aT yt に関して， TDMI を最大化することを考える．この時，式 (27) は次. 3.
(4) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. のようになる．. ここで, 行列 Kx を第 (n, m) 要素が,. Mx(k) st = t−1. aT Σyt yt a 1 . log2 T 2 a Σyt yt |x(k) a. aT Σyt yt a = 1 という拘束条件のもとで式 (28) を最大化す ることを考える．ラグランジュの未定乗数 λ を用いること で以下の式が得られる．. (33). で定められる N × N の対称行列として定義する. この行列 はグラム行列と呼ばれる. グラム行列を用いると, データ ˆ xx は xn の再生核ヒルベルト空間における共分散作用素 Σ 次の式で表される.. a Σyt x(k) Σ−1 (k) Σ (k) (k) xt−1 xt−1 xt−1 yt t−1. = λΣyt yt a.. (29). これも Dy 個の固有値とこれに対応する固有ベクトル Dy {λi , ai }i=1 ，(λ1. Kx,nm = kx (xn , xm ). (28). t−1. ˆ xx = K ˜ 2, Σ x. (34). ˜ x = Kx − 1N Kx − Kx 1N + 1N Kx 1N である． ただし K. ≥ · · · ≥ λDy ≥ 0) を持つ一般化固有値問題. 1N は全ての要素が 1/N という値を取る N × N 行列であ. となる．この問題は正準相関分析（Canonical Correlation. る. データ yn についても正定値カーネル関数 ky (yn , ym ). (k). Analysis，以下 CCA と略す）と等価である．xt−1 の線形. によって写像した再生核ヒルベルト空間を用いて考える. 行. 写像についても考えると次のようになる．. 列 Ky を第 (n, m) 要素が,. Σx(k) yt Σ−1 yt yt Σyt x(k) b = λΣx(k) x(k) b. t−1. t−1. t−1. (30). Ky,nm = ky (yn , ym ). t−1. (35). 式 (29) と (30) は共通の固有値を持つ．これらの固有値を. ˆ yy および Σ ˆ xy は同 の対称行列とすると, 共分散作用素 Σ. 用いることで，以下の式を得る．. 様に次の式で表される.. ∑ Dy. Mx(k) yt =. M. t−1. i=1. ここで M. {i} (k). xt−1 yt. ∑1 Dy. {i} (k). xt−1 yt. =. i=1. 2. log2. 1 , 1 − λi. (31). は i 番目の部分空間の TDMI の値である．. このようにガウス分布の仮定のもとでは，TDMI は CCA の固有値の評価と同等となる．. ˆ yy = K ˜ 2, Σ y. (36). ˆ xy = K ˜ xK ˜y, Σ. (37). ˜ y = Ky − 1N Ky − Ky 1N + 1N Ky 1N である． ただし K kx (xn , xm ) と ky (yn , ym ) は異なる正定値カーネル関数で あってもかまわない. このとき, xn が与えられたときの yn. この時点で，TE と TDMI を一般化固有値問題という同. の条件付き共分散作用素は,. じ視点から評価できるようになる．TDMI と TE はほぼ同 (l). じような構造を持つことが分かるが，yt−1 を事前情報とし て用いるかが異なる点となる．より詳しく言うと，確率密 (k). ˆ yy|x = Σ ˆ yy − Σ ˆ yx Σ ˆ −1 Σ ˆ Σ (38) xx xy ( 2 ) ( 2 )−1 ˜ xK ˜y ˜ + ϵIN − K ˜yK ˜x K ˜ + ϵIN K = K y x. 度分布としてガウス分布を仮定すると，TDMI は xt−1 と (l) yt の間のクロス相関を yt−1 を考慮せずに測っていること (l) (k) になる．一方，TE は yt−1 が与えられた下での xt−1 と yt (l) との偏相関を測っていることになる．つまり，yt−1 の影響. が大きければ，TE と TDMI との差が大きくなる．. (39) ˆ xx や Σ ˆ yy は逆作用素を考える必要がある場合が となる. Σ ˆ xx あるが, 一般にこれらの作用素は非可逆である. そこで Σ ˆ yy に対して正則化を行う. やΣ カーネル PCCA は次の一般化固有値問題で定義される.. 4. カーネル化. ˆ (k) (l) Σ ˆ −1 Σ (k) yt x |y. これまで議論してきたように，多次元時系列の因果指標. t−1. t−1. (k). (l). xt−1 xt−1 |yt−1. ˆ (k) Σ (l) α x yt |y t−1. t−1. ˆ = λΣ (l) α, yt yt |y. の算出は CCA の一種である PCCA に帰着される．ここで. (40). t−1. はカーネル PCCA を用いた因果指標を提案する. 主成分分析（PCA）や CCA にカーネルトリックを適用 したカーネル PCA やカーネル CCA では，正定値カーネ ル関数 kx (xn , xm ) による N 個のデータ xn の写像. Φ : xn → kx (x, xn ) の線形結合作用素. ∑N n=1. (32). αn kx (x, xn ) で構成されるベク. トル空間, 再生核ヒルベルト空間において通常の PCA や. CCA を行う. カーネル PCCA も同様に, 再生核ヒルベルト. ただし. ( 2 ) ˜ + ϵIN ˆ K Σ (l) = yt yt yt |yt−1 ( )−1 ˜y K ˜ (l) K ˜ 2(l) + ϵIN ˜ −K K t. yt−1. ˆ (k) (k) (l) = Σ x x |y t−1. t−1. t−1. ˜ (k) K ˜ (l) −K x y t−1. t−1. yt−1. ( ) ˜ 2 (k) + ϵIN K x. (l). yt−1. ˜y , K t. (41). t−1. ( )−1 2 ˜ ˜ (l) K ˜ (k) , (42) Ky(l) + ϵIN K y x t−1. t−1. t−1. 空間での PCCA として定義される.. c 2014 Information Processing Society of Japan ⃝. 4.
(5) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. ˆ (k) (l) = K ˜y K ˜ (k) Σ t yt xt−1 |yt−1 xt−1 ( )−1 ˜y K ˜ 2(l) + ϵIN ˜ (l) K ˜ (k) , (43) ˜ (l) K K −K t y y x y t−1. t−1. t−1. t−1. ˆ (k) ˜ ˜ Σ (l) = K (k) Kyt xt−1 yt |yt−1 xt−1 ( )−1 2 ˜ ˜ ˜ y . (44) ˜ ˜ (l) K − Kx(k) Ky(l) Ky(l) + ϵIN K t y t−1. t−1. (l). (k). ′. T. t−1. t−1. yt , yt−1 , xt−1 の 全 て に 対 し て 線 形 な カ ー ネ ル 関 数 ′. k (x, x ) = x x を用いた場合は, 通常のデータ空間で の PCCA と等価である. この一般化固有値問題を解いて得 られる固有値 λi を用いて, 変数 x の部分空間から y の部 分空間への非線形な TE は次のようになる． {i} Tx→y =. 1 ||X′ − ZA||2 + ηz AT A, N. ここで ηz は正則化のパラメータである．この様な，係数 の 2 乗に比例する正則化項を付加した回帰分析のことを， リッジ回帰と呼ぶ．. ˇ を求める．A に関す X′ = [X Y] として A の推定値 A る J の導関数を 0 とすると，. ˇ = (Σzz + ηz I)−1 [Σzx Σzy ] A =. ˜ −1 Σ zz. ∑. [Σzx Σzy ] .. (49) (50). ˜ zz = Σzz + ηz I と置いた．得られた ただし，式 (50) で，Σ ˇ による X，Y の予測値と実際の値との間の誤差の分散 A. (45) {i} Tx→y. の総和は 2. 偏共分散行列は以下のようになる．. [. Σ{xy}{xy}|z =. {i} Tx→y .. (46). i. 5. 正則化 多次元データに対して回帰分析や CCA を適用する時の 問題点として，データの各次元の間の相関が強い場合や，. (データの次元数)>(サンプル数) である場合に，計算に使. Σxx. Σxy. ]. Σyx Σyy [ ] Σxz ˇ +A ˇ T Σzz A ˇ −2 A Σyz [ ] Σxx|˜z Σxy|˜z = . Σyx|˜z Σyy|˜z. つの変数 x 全体から y 全体への非線形な TE となる.. Tx→y =. (48). 共分散行列の事を偏共分散行列と呼び，正則化した場合の. 1 1 log2 . 2 1 − λi. このとき非ゼロの固有値から算出された. J=. (51). (52). ここで，. ˜ −1 ˜ −1 ˜ −1 Σxx|˜z = Σxx − 2Σxz Σ zz Σzx + Σxz Σzz Σzz Σzz Σzx ,. う分散共分散行列がランク落ち，あるいはそれに近い状態. ˜ −1 Σzy + Σxz Σ ˜ −1 Σzz Σ ˜ −1 Σzy , Σxy|˜z = Σxy − 2Σxz Σ zz zz zz. になり逆行列計算が不可能または不安定になることがあげ. ˜ −1 Σzx + Σyz Σ ˜ −1 Σzz Σ ˜ −1 Σzx , Σyx|˜z = Σyx − 2Σyz Σ zz zz zz. られる．そこで，多次元 TE に対し正則化を行う [7]．本章. ˜ −1 Σzy + Σyz Σ ˜ −1 Σzz Σ ˜ −1 Σzy . Σyy|˜z = Σyy − 2Σyz Σ zz zz zz. では. (m) Xt. を，ある時刻 ts から te までの埋め込みベクト (1). ルを並べた行列とする．ただし，Xt = Xt. である．. ここまで，リッジ回帰を利用することにより，変量 Z に おける共線性などによる分散共分散行列 Σzz の逆行列計. (m). Xt. = ( x(m) ts. (m). xts +1. ···. (m). xt e. )T .. 算の不安定性を緩和し，その場合の偏共分散行列を導出し た．しかし，変量 X，Y に関する正則化はできていない. 5.1 偏正準相関分析の正則化. ため，Σxx ，Σyy がランク落ちしていた場合，PCCA を解. 多次元 TE（あるいは GC）は PCCA で表現されるため，. くことはまだ出来ない．この様な場合にも使える CCA の. 多次元の TE の正則化を行うには，PCCA に対して正則化. 正則化手法として，Canonical Ridge という手法がある．. を行えば良いことになる．. Canonical Ridge は以下の様な，相関係数を最大化する射. まずは，正則化した場合の偏共分散行列について考える．. X ′ ，Z をそれぞれ N × Dx′ ，N × Dz 行列とし，係数行列 A を Dz × Dx′ 行列，残差 ϵx′ を N × Dx′ 行列とした場合 の線形回帰式は，以下のようになる．. X′ = ZA + ϵx .. 影を見つける問題として表現される．. ρ =√. aT Σxy b √ , aT (Σxx + ηx I) a bT (Σyy + ηy I)b. (53). ηx ，ηy は正則化のパラメータ (ηx , ηy ≥ 0) である． (47). 次に正則化を行った回帰の場合の偏共分散行列を求める． 正則化とは誤差関数に正則化項を付加することにより，係 数が必要以上に大きくなる事などを防ぐ効果を得るもので. 偏共分散行列を，正則化を行った場合のものに変更した. PCCA に対して，式 (53) と同じように記述すると以下の 様になる．. ρ =√. ある．また，この場合は Z の分散共分散行列の逆行列計算. aT Σxy|˜z b √ . aT (Σxx|˜z + ηx I)a bT (Σyy|˜z + ηy I)b. (54). を安定化する効果がある．正則化項を付加した誤差関数は. ラグランジュの乗数 λ，ν を用いて，次の最大化問題を考. 以下のようになる．. える．. c 2014 Information Processing Society of Japan ⃝. 5.
(6) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. √ ) λ( T J = a Σxy|˜z b − a (Σxx|˜z + ηx I)a − 1 2 √ ) ν( T − b (Σyy|˜z + ηy I)b − 1 . 2. Σx(m) yt |˜y(n) = Σx(m) yt. T. t−1. T. T. となり，これらの式に左からそれぞれ a ，b を掛けて式 √ √ (55) に代入すると， λ = ν = aT Σxy|˜z b が成り立つ．a，. b のいずれかを消去すると，次の一般化固有値問題を得る．. t−1. − Σx(m) y(n) (Σy(n) y(n) + ηy(n) I)−1 Σy(n) yt , t−1 t−1. (55). a，b それぞれに関する J の導関数を 0 とすると， √ ∂J = Σxy|˜z b − λ(Σxx|˜z + ηx I)a, (56) ∂a √ ∂J = Σyx|˜z a − ν(Σyy|˜z + ηy I)b, (57) ∂b. t−1. t−1 t−1. t−1. (66). t−1. Σyt x(m) |˜y(n) = Σyt x(m) t−1. t−1. t−1. − Σyt y(n) (Σy(n) y(n) + ηy(n) I)−1 Σy(n) x(m) . t−1 t−1. t−1. t−1. t−1. (67). t−1. ηx(m) ，ηy(n) ，ηyt は そ れ ぞ れ 正 則 化 の パ ラ メ ー タ で あ t−1. t−1. る（ηx(m) , ηy(n) , ηyt ≥ 0）．式 (60)，あるいは式 (61) の t−1. t−1. 一般化固有値問題から固有値 λ1 ≥ λ2 ≥ · · · ≥ λM ≥. 0 (M = min{m × Dx , Dy }) が得られる． 第 i 正準空間同士の TE，x 全体から y 全体への TE は それぞれ式 (45)，式 (46) と同様に表現できる．. ˜ xx|˜z a, ˜ −1 Σyx|˜z a = λΣ Σxy|˜z Σ yy|˜ z. (58). ˜ −1 Σxy|˜z b = λΣ ˜ yy|˜z b. Σyx|˜z Σ xx|˜ z. (59). 6. 確率的解釈とベイズモデル この章では因果指標計算に用いる PCCA と等価な確率. ˜ xx|˜z = Σxx|˜z + ηx I, Σ ˜ yy|˜z = Σyy|˜z + ηy I． ただし，Σ. 的生成モデルを与える．これにより確率的観点からのモデ ル拡張が可能になる．その拡張としてモデルパラメタに事. 5.2 多次元 Transfer Entropy の正則化. 前分布を仮定しパラメタを分布推定するベイズ推定のモデ. ここからは，先程導出した正則化 PCCA を，実際に多. ルを提案する [8]．. 次元の TE に用いた時の定式化を行う．前節まで扱ってき た PCCA は，Z の情報を取り除いた状態で，X と Y の相 関を最大化するという問題であった．データにガウス分布 を仮定したときの多次元 TE では，変量 X から Y への因 (m). (n). 6.1 生成モデル 影響を除きたい第三変数 xn ，潜在変数 zn から同時に線 形回帰する以下の式に表される生成モデルを考える．. 果を見るときは，Yt−1 の情報を取り除いた状態で，Xt−1. zn ∼ N (0, Idz ),. と Yt の相関を最大化する問題を解くことになる． 以上のことから正則化を行った PCCA に対して， (n). (m). X → Xt−1 , Y → Yt , Z → Yt−1 ,. t−1. (m). (n). xt−1 xt−1 |˜ yt−1. zn ∼ N (0, Idz ),. Σx(m) yt |˜y(n) wyt t−1. t−1. ˜ = λΣ (n) wyt , y yt yt |˜. (60). t−1. ˜ −1 Σx(m) yt |˜y(n) Σ t−1. t−1. (n) yt yt |˜ yt−1. t−1. t−1. t−1. t−1. (61). ( C= (62). ˜ Σ (n) = Σ (n) + ηyt I, yt yt |˜ y yt yt |˜ y. (63). t−1. t−1. t−1. t−1. t−1. t−1. t−1. t−1. t−1. t−1. n=1. t−1 t−1. t−1. t−1. t−1. t−1. − Σyt y(n) (Σy(n) y(n) + ηy(n) I)−1 Σy(n) yt , t−1 t−1. +. Wz1. ) (. Wz2. T. Wz1 Wz2. T. ). ((. µ1 µ2. ) ( −. yn1. ) (. yn2. +. Wx1 Wx2. ). ) xn. t−1. Σyt yt |˜y(n) = Σyt yt t−1. ) (. 0 Ψ2. ∑ ∂L =− C −1 ∂µ N. − Σx(m) y(n) (Σy(n) y(n) + ηy(n) I)−1 Σy(n) x(m) , (64) t−1 t−1. Ψ1 0. と置くと，. Σx(m) x(m) |˜y(n) = Σx(m) x(m) t−1. の y の共分散を. t−1. ˜ (m) (m) (n) = Σ (m) (m) (n) + η (m) I, Σ x x |˜ y x x |˜ y x t−1. (69). 今，対数尤度を L と置き，潜在変数 z を積分消去した時. t−1. ただし，次の様に変数を置き換えた．. t−1. ynm ∼ N (Wzm zn + µm , Ψm ) . に帰着されることを示す．. Σyt x(m) |˜y(n) wx(m). ˜ (m) (m) (n) w (m) . = λΣ x x |˜ y x t−1. PCCA を用いて計算出来ることを示す．そのために，提案 モデルが確率的正準相関分析の式. 一般化固有値問題で表現される．. t−1. (68). このモデルの最尤推定解 arg max log p(y|x; Wx , Wz , Ψ) が. と置き換えると，多次元 TE を正則化したものは，以下の. ˜ −1 Σyt x(m) |˜y(n) Σ (m). ynm ∼ N (Wxm xn + Wzm zn + µm , Ψm ) .. t−1. (65). が成り立つ．C は正定値なので，対数尤度を最大にする µ は偏微分が 0 になる点である．従って. µm = y m − Wxm x. (70). t−1. となる．各サンプルからサンプル平均を引いたものを. c 2014 Information Processing Society of Japan ⃝. 6.
(7) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. y en1 = yn1 − y 1 のように書くと，式 (70) を代入し， ) (( ) ( ) N ∑ y en1 Wx1 ∂L −1 T T = C x en − x en x en ∂Wx y en2 Wx2. であれば推論に影響しない．ここで Gamma(a, b) はガ. となる．この式はサンプルが入力空間を張る時 L が Wx の. くなるように小さいものがこのましい．しかしながら逆. 負定値二次形式になることが示せるので，偏微分が 0 とな. Wishart 分布の定義から ν0m > dm − 1 となる．ARD 事前. る Wx で L は最大化される．従って，. 分布により重要度の低い列が 0 に向かうため，十分大きな. m ンマ分布， IW(ν, ( ) K) は逆 Wishart 分布である．W =. m は W m の k 列目を表す．ハ Wxm Wzm であり，W:,k イパーパラメタ a0 ，b0 ，ν0m ，K0m は事前分布のすそが広. n=1. Wxm = Σmx Σ−1 xx. dz をとれば適切な潜在変数 z の次元が推定できる．. (71). 推論は変分ベイズ法により行う．近似した事後分布を. と な る ．こ の 式 を 生 成 モ デ ル の 式 (68) に 代 入 す る と, 確 率 的 正 準 相 関 分 析 の 式 (69) に 入 力 と し て. m y′ n. =. q(Z, Θ) = q(Z). y enm − Σmx Σ−1 en を代入したものと等価になる．これらの xx x. 2 ∏. (. m. m. q(Ψ )q(α ). m=1. dm ∏. ). q(wjm ). (75). j=1. と置き，各分布を逐次的に更新する．ここで wjm は W m. 共分散は. の j 番めの行である．各事後分布は. N 1 ∑ ′ m1 ′ m2 T yn yn = Σm1 m2 − Σm1 x Σ−1 xx Σxm2 N. q(zn ) = N (µzn , Σzn ),. n=1. = Σm1 m2 |x. (72). と，偏共分散になるため, 最尤解は確率的正準相関分析の 最尤解で共分散を偏共分散で置き換えたもの，すなわち. q(Ψm ) = IW(νm , Km ), q(wjm ) = N (µwjm , Σwjm ), ∏ q(αm ) = Gamma(am , bmk ). PCCA に帰着されることが示された．式で表現すると，最 尤解は以下に表される．. の形になり，パラメタの更新式は. (. Wxm = Σmx Σ−1 xx ,. Σzn = I +. Wzm = Σmm|x Udmz Mm ,. ∑. )−1 ⟨(Wzm )T (Ψm )−1 Wzm ⟩. ,. m. Ψm = Σmm|x − Wzm Wzm T , µm = y m − Wxm x,. (76). k. µzn =Σzn (73). ∑( ⟨(Wzm )T ⟩⟨(Ψm )−1 ⟩ynm m. ) − ⟨(Wzm )T (Ψm )−1 Wxm ⟩xn ,. ここで Udmz は d 列に d 番めの固有ベクトルを並べた行列，. Pd は対応する固有値を d 個対角に並べた対角行列，Mm は. M1 M2T. = Pdz を満たしスペクトルノルムが 1 より小さ. い任意の行列である． 以下の議論では簡潔性のためデータは平均が 0 になって いると仮定し，µ の推論は行わない．. 6.2 ベイズモデル この節では，Wang [9] の手法にしたがい，生成モデルの パラメタ Ψ, W に事前分布を仮定し，パラメタの事後分布 をベイズ推定するモデルを考える．これによりロバストな 推定，モデル選択が行われる．W に列毎に ARD 事前分 布 [10] を仮定し，Ψ に逆 Wishart 分布を仮定する．生成モ デルは以下の式に表される．. αkm. ( m. Σwjm = diag⟨α. ⟩+⟨(Ψm )−1 j,j ⟩. m µwjm =⟨(Ψm )−1 j,: ⟩Y. ∑. (. XT. m ⟨(Ψm )−1 j,l ⟩⟨Wl,: ⟩. m bmk =b0 + ⟨||W:,k ||⟩/2,. +. Wzm zn , Ψm ),. (74). ここで第三変数の事前分布 p(x) は各サンプルで p(x) > 0. c 2014 Information Processing Society of Japan ⃝. XX T X⟨Z⟩. T. ))−1 ,. ⟨Z⟩X T ⟨ZZ T ⟩ ). ZT (. T. XX T X⟨Z⟩. ⟨Z⟩X T ⟨ZZ T ⟩. ) ,. am =a0 + dm /2,. zn ∼ N (0, Idz ), ∼. (. l̸=j. Ψm ∼ IW(ν0m , K0m ), N (Wxm xn. νm =ν0m + N,. −. ∼ Gamma(a0 , b0 ),. m W:,k ∼ N (0, (αkm )−1 Idm ),. ynm. Km = K0m + Y m (Y m )T ) ( XX T XZ T m (W m )T ⟩ + ⟨W ZX T ZZ T ( ) − Y m X T ⟨Z T ⟩ ⟨(W m )T ⟩ ( ) X m − ⟨W ⟩ Y m, ⟨Z⟩. (77). の形になる．ここで，diag⟨αm ⟩ は k 番めの対角要素が. ⟨αm k ⟩ となる対角行列である． 7.
(8) Vol.2014-CG-157 No.16 Vol.2014-CVIM-194 No.16 2014/11/21. 情報処理学会研究報告 IPSJ SIG Technical Report. 6.3 近似されたベイズモデル. aαm =a0 + dm /2,. 前の節で提案したベイズモデルはノイズ精度行列 Ψ の. b αm =b0 + ⟨(W m )T W m ⟩k,k /2, k. 推論に大きな計算量が必要であり，またサンプル数が小さ. aτ m =a0 + N dm /2, ( ( 1 bτ m =b0 + Tr Y m (Y m )T 2 ( ) ) − 2Y m X T ⟨Z T ⟩ ⟨(W m )T ⟩ ( ( ))) T XX T X⟨Z⟩ m T m + Tr ⟨(W ) W ⟩ . (82) ⟨Z⟩X T ⟨ZZ T ⟩. いときに逆 Wishart 事前分布の影響が大きくなる．そこで. Klami ら [11] の手法に従って，ノイズを等方ノイズと非共 有の潜在変数からの線形回帰で近似する以下の生成モデル を考える．. zn ∼ N (0, Idz ), znm ∼ N (0, Idzm ), ( ) ynm ∼ N Wxmxn+Amzn +B mznm , (τ m )−1Idm .. (78). 7. おわりに. 非共有潜在変数 znm を積分消去すると，前の節の生成モ. 本論文では，まず，Transfer Entropy，Continuous Trans-. デル式 (74) において Ψm = B m (B m )T + (τ m )−1 Idm とし. fer Entropy，Granger Causality といった従来の因果指標. たものと等価になる. 従ってこのモデルは共分散を低ラ. を紹介し，これらの関係性を明らかにした．特に時系列確. ンク近似したものとみなせる. A(と B を同時に推論する ) A(1) B (1) 0 ため，これらをまとめて Wz = , A(1) 0 B (1). 率変数がガウス分布に従うときに，Transfer Entropy が偏. (. W = αkm m W:,k. ). Wx. とかき，以下のモデルを考える．. Wz. をカーネル化することで，様々な計量を考慮可能な因果指 標を提案した．また，偏正準相関分析の正則化により，少. ∼ Gamma(a0 , b0 ), ∼. 正準相関分析に帰着されることを示した．偏正準相関分析. サンプルでも安定して計算可能な手法を提案した．さら に，偏正準相関分析を確率的に解釈をし，確率的偏正準相. N (0, (αkm )−1 Idm ),. 関分析という新たなモデルを提案をした．これにより確率. τ m ∼ Gamma(a0 , b0 ),. 的観点からのパラメタ推定やモデル拡張が可能になった．. zn ∼ N (0, I(dz +dz1 +dz2 ) ), ( ) ynm ∼ N Wxm xn + Wzm zn , (τ m )−1 Idm .. (79). この表現によりモデルパラメタ数が減り, さらなるロバス ト性が期待される．このモデルにおいてもハイパーパラメ. 参考文献 [1] [2]. タ a0 ，b0 は小さいものが好ましい． 近似された事後分布を. q(Z, Θ) = q(Z). ∏. [3]. (q (τ m ) q (αm ) q (W m )). (80). m. [4]. と置くと，変分ベイズ推論により各々の事後分布は. q(Z) =. ∏. N (µzn , Σz ),. n. q(W m ) =. ∏. [5]. m , ΣW m ), N (µWd,:. d. m. q(α ) =. ∏. [6]. Gamma(aαm , bαm ), k. k. q(τ m ) = Gamma(aτ m , bτ m ). (81). の形になり，パラメタ更新則は以下の通りとなる．. (. (. ΣW m = diag⟨α ⟩+⟨τ ⟩ m. µW m =Y m ( Σz =. ( X. I+. ∑. ( ∑. ⟨Z⟩=Σz. T. ⟨τ. T. XX T X⟨Z⟩. m. ⟨Z ⟩ T. m. ). ))−1. ⟨Z⟩X T ⟨ZZ T ⟩. [7] [8]. , [9]. ,. )−1. ⟩⟨(Wzm )T Wzm ⟩. [10]. ,. m. ( ⟨τ m ⟩ ⟨(Wzm )T ⟩Y m−⟨(Wzm )T Wxm ⟩X. ) ). [11]. ,. Schreiber, T.: Measuring Information Transfer, Phys. Rev. Lett., Vol. 85, No. 2, pp. 461–464, (2000). Kaiser, A. and Schreiber, T.: Information transfer in continuous processes, Physica D Nonlinear Phenomena, Vol. 166, pp. 43–62 (2002). Granger, C. W. J.: Investigating Causal Relations by Econometric Models and Cross-Spectral Methods, Econometrica, Vol. 37, No. 3, pp. 424–438 (1969). Shibuya, T., Harada, T. and Kuniyoshi, Y.: Causality Quantification and Its Applications: Structuring and Modeling of Multivariate Time Series, KDD, pp. 787–796 (2009). Shibuya, T., Harada, T. and Kuniyoshi, Y.: Reliable index for measuring information flow, Phys. Rev. E, Vol. 84, p. 061109 (2011). Rao, B. R.: Partial canonical correlations, Trabajos de estadistica y de investigacion operativa, Vol. 20, No. 2-3, pp. 211–219 (1969). Yamashita, Y., Harada, T. and Kuniyoshi, Y.: Causal Flow, IEEE TMM, Vol. 14, No. 3, pp. 619–629 (2012). Mukuta, Y. and Harada, T.: Probabilistic Partial Canonical Correlation Analysis, ICML, pp. 1449–1457 (2014). Wang, C.: Variational Bayesian Approach to Canonical Correlation Analysis, Trans. Neur. Netw., Vol. 18, No. 3, pp. 905–910 (2007). Neal, R. M.: Bayesian Learning for Neural Networks, PhD Thesis, the University of Toronto (1995). Klami, A., Virtanen, S. and Kaski, S.: Bayesian Canonical Correlation Analysis, J. Mach. Learn. Res., Vol. 14, No. 1, pp. 965–1003 (2013).. m. c 2014 Information Processing Society of Japan ⃝. 8.
(9)
関連したドキュメント
ユーザー情報のダウンロード エラー内容 要因① ウイルスソフト関連 要因② Proxyサー バー環境. 要因③
We found that the use of algorithms for causality assessment, the grades and expressions used to describe causality, and criteria for determining whether reactions were ADRs or
最も偏相関が高い要因は年齢である。生活の 中で健康を大切とする意識は、 3 0 歳代までは強 くないが、 40 歳代になると強まり始め、
[r]
目標 目標/ 目標 目標 / / /指標（ 指標（ 指標（KPI 指標（ KPI KPI KPI）、実施スケジュール ）、実施スケジュール ）、実施スケジュール ）、実施スケジュールの の の の設定
添付資料 4 SDC 3/INF.10: Information collected by the intersessional Correspondence Group on Intact Stability regarding second generation intact
この標準設計基準に定めのない場合は，技術基準その他の関係法令等に
様々な国の子供の死亡原因とそれに対する介入・サービスの効果を分析すると、ミレニ アム開発目標 4
