混合確率的偏正準相関分析
Mixture of Partial Canonical Correlation Analysis
河野圭祐
横山裕樹
森裕紀
浅田稔
大阪大学大学院工学研究科 知能・機能創成工学専攻
The Japanese Society for Artificial IntelligenceIn this article, we propose Mixture of Probabilistic Partial Canonical Correlation Analysis (MPPCCA) to analyze data that include multiple causal relationships. Partial Canonical Correlation Analysis (PCCA) is a method for multivariate statistics. PCCA evaluates a correlation between two multivatiate while considering third multivariate, and has been used as a calculation method of Granger Causality which is one of the causal indexes.
However, PCCA is not applicable to analyze data which contaion some different partial canonical correlations. In this paper, we addressed the problem by proposing a Mixture of Probabilistic Partial Canonical Correlation (MPPCCA) and deriving a parameters estimation method by EM algorithmn.
We show an experiment with artificial data which indicates our mehod calculate more accurately than exiting method. An actual data is also analyzed by the proposed mothod. The data was divided by into partial correaltion.
1.
はじめに
偏正準相関分析[1]は,ある多変量の影響を除いた上で,2 つの多変量の相関が大きくなる低次元空間への射影を求める多 変量解析手法であり,偏相関分析の多変量への拡張といえる. 偏正準相関分析は,時系列データの因果解析などに応用されて いる[2].また,高次元少サンプルなデータに対しても安定し て計算を行うことを目的とした,確率的な解釈も与えられてい る[3]. 一方で,データ中の偏正準相関関係は全て一定であるとは 限らず,複数の関係が混在している可能性がある.例えば,時 系列データの因果解析では,多変量間の因果関係が時間ととも に切り替わる状況(図1)が考えられる.偏正準相関分析では, このようなデータを正しく解析することはできない. そこで,本研究では確率的偏正準相関分析[3]の混合モデル を提案し,EMアルゴリズム[4]による,モデルパラメータの 更新式を導出した.提案手法では,複数の偏正準相関関係が混 在しているデータについて,同一の偏正準相関関係であること を基準としたデータ分割をすることができる.提案手法の検証 として,人工的に発生させた,複数の偏正準相関関係が混在し たデータを用いた実験によって,提案手法と既存手法[5]を比 較し,提案手法が既存手法と比べて,データを発生させるとき に用いた偏正準相関関係を正しく推定できることを確認した2.
関連研究
本章では本研究に関連する多変量解析手法について説明す る.さらに,偏正準相関分析とグレンジャー因果[6]との関係 についても簡単に述べる.2.1
正準相関分析
正準相関分析[7, 8]は二つの多変量の相関関係を分析する際 に用いられる手法であり,脳波の解析や天気の予報など幅広く 応用されている.正準相関分析では二つの多変量の相関が最大 になる低次元空間を求め,それぞれの多変量をその空間へ射影 し,相関係数を計算することで,二つの多変量間の相関関係を定 量化する.二つの多変量y(1)= (y(1)1 , y2(1), ..., ym(1)1) T ∈ Rm1, 連絡先:河野圭祐,[email protected] !! !! !! !! !! !!"#"#!
図1: パターンごとに異なる関係をもつ時系列 y(2) = (y(2) 1 , y (2) 2 , ..., y (2) m2) T ∈ Rm2 の間の正準相関分析は以 下の一般化固有値問題に帰着する[9]. (Σ12TΣ11−1Σ12− ρ2Σ22)u(2) = 0 (Σ21TΣ22−1Σ21− ρ2Σ11)u(1) = 0 ここでΣab ∈ Rda×dbはサンプルy(a)∈ Rda,y(b)∈ Rdb の共分散行列である.また,ρは正準相関係数,つまり相関関 係の強さである.各固有ベクトルが低次元空間への射影と対応 しており,対応する固有値が射影した空間での相関,つまり正 準相関に対応する.2.2
偏正準相関分析
偏正準相関分析[1]はある一つの多変数の影響を除いた上で の二つの多変量の正準相関を求める手法であり,条件を統制し た上での正準相関分析や後述するグレンジャー因果[6]の計算 [2]などに用いられている.多変量x = (x1, x2, ..., xd)T の影 響を除いた上で,二つの多変量y(1) = (y1(1), y2(1), ..., ym(1)1) T , y(2)= (y(2)1 , y2(2), ..., y(2)m2) T の間の変正準相関分析は以下の一 般化固有値問題に帰着する[9]. (ΣT12|xΣ−111|xΣ12|x− ρ 2 Σ22|x)u (2) = 0 (ΣT21|xΣ−122|xΣ21|x− ρ2Σ11|x)u(1) = 0 (1) ここで,Σab|x = Σax− ΣaxΣ−1xxΣxbである.またρは偏正 準相関係数であり,xの影響を除いた上での,二つの多変量 y(1),y(2)の相関関係の強さを表す.1
The 29th Annual Conference of the Japanese Society for Artificial Intelligence, 2015
図2: 確率的偏正準相関分析のグラフィカルモデル[3]
2.3
因果指標と偏正準相関分析
偏正準相関分析を用いることで,因果指標の一つであるグレ ンジャー因果[6]を求めることができる[2, 10].グレンジャー 因果は時系列データ間の因果関係を示す量であり,経済学や脳 神経科学などの分野におけるデータ解析の手法として広く用い られている.グレンジャー因果は時系列データの予測精度に基 づく,統計的な因果指標であり,実際の因果関係とは必ずしも 一致しない. 二つの多次元時系列データX,Y が定常かつ平均0のとき, Y からXのグレンジャー因果は,以下の一般化固有値問題を 解くことで求まる[9]. ( ΣT XtYt−k(m)|Xt−k(m)Σ −1 XtXt|Xt−k(m) Σ XtYt−k(m)|X(m) t−k −ρ2 Σ Yt(m)−kYt(m)−k|Xt(m)−k ) a = 0 (2) ここで,Σab|c = Σac− ΣacΣcc−1Σcbである.また,Xt(m)−1 = (XtT−1XtT−2...XtT−m)T,Yt(m)−1 = (YtT−1YtT−2...YtT−m)Tは埋め 込みベクトルである.ここで,式(2)は,偏正準相関関係を求 めたい二つの変数をXt,Yt(m)−1 とし,影響を除きたい変数を Xt(m)−1とした場合の偏正準相関分析に等しい. このときグレンジャー因果はそれぞれの固有値ρiを用いて 以下のように求まる. GY→X = 1 2log2 1 1− ρ2 i (3) 大きい固有値に対応したグレンジャー因果の値は大きくなる. 本研究ではこれらのグレンジャー因果のうち,最大のものをグ レンジャー因果の値として,解析を行う.2.4
確率的偏正準相関分析
Mukuta et al.[3]は,図2に示すグラフィカルモデルにより 偏正準相関分析に確率的な解釈を与えた.偏正準相関を考えた い多変量をy(1)∈ Rd1,y(2)∈ Rd2,影響を除きたい多変量を x∈ Rdxとするとき,各確率変数の確率密度関数は以下のよ うになる. p(tn) = N (0, Id) p(yn(1)|tn) = N (Wx(1)xn+ W (1) t tn+ µ(1), Ψ(1)) p(yn(2)|tn) = N (Wx(2)xn+ Wt(2)tn+ µ(2), Ψ(2)) 観測変数y(1),y(2)は共通因子tおよび,xからも影響を受 けており,Wt, Wxはそれらの強さを表す.偏正準相関分析で 求めたdt個の固有値を対角成分にもつ対角行列をPd,対応 する固有ベクトルを並べた行列をU1d,U2dとしたとき,パラ メータの最尤解は以下のようになる. Wx(m) = ΣmxΣ−1xx Wt(m) = Σmm|xUd(m)t Mm Ψ(m) = Σmm|x− Wt(m)W (m)T t µ(m) = y(m)− Wx(m)x ただし,m = 1, 2であり,Σmm|x= Σmm− ΣmxΣ−1xxΣxmで ある.またy(m)= 1 N ∑N n=1y (m) n ,x =N1 ∑N n=1xnを表す.2.5
空間分割による因果パターン抽出
朴ら[5]は偏正準相関分析を用いた因果性評価手法において, 因果性の原因側の過去情報と結果側の過去,現在情報で構成さ れる特徴空間を想定し,その空間におけるデータ点をクラス タリングすることで,因果性の高いパターンを抽出した.朴ら [5]の手法では,埋め込みベクトルXt(m)−1,Yt(m)−1 を用いて特徴 点Ψt= (Xt, Xt(m)−1, Yt(m)−1)が構成される.この特徴点Ψには 因果性の原因側の過去情報と結果側の過去,現在情報が含ま れている.朴ら[5]はこれらの特徴点の集合で構成される空間 において,ユークリッド距離が近い点同士は同じ因果関係で説 明できると仮定し,この空間上でK-meansクラスタリングを 行い,そのクラスタごとにグレンジャー因果を計算することに よって因果性を評価した.しかし,この手法はユークリッド距 離が近いデータ点同士は同じ因果関係で説明できると仮定して いるため,特徴量は似ているが全く異なる因果関係を持ってい るデータに対して正しく因果性を評価することができない.3.
確率的偏正準相関分析の混合モデル
本章では,確率的偏正準相関分析[3]の混合モデルを提案し, EMアルゴリズムによるパラメータおよび潜在変数の最尤推 定を行う.提案手法では,同一の偏正準相関関係または同一の グレンジャー因果で説明可能な部分データ集合を抽出し,その データ集合ごとに偏正準相関の評価をする.3.1
生成モデル
偏正準相関関係を求めたい二つの多変量をy(1) ∈ Rm1, y(2)∈ Rm2とし,影響を取り除きたい多変量をx∈ Rdxとす る.また,潜在変数としてtn∈ Rdtおよびznk∈ {0, 1}を導 入する.ただしmin(d1, d2)≥ dtである.tnはy(1),y(2)間 の共通因子を表す.また,znkはn番目のデータ点がどのクラ スタから生成されたものかを表しており,観測データからznk を求めることで,各点がどのクラスタから生成されたものかを 推定することができる.∑kznk = 1であり,zkはどれか一 つのkに対してzk= 1で,それ以外のkに対してzk= 0で ある. p(tn) =N (tn| 0, Idt) p(y(1)n | tn, xn; W (1) xk, W (1) tk , µ (1) k , Ψ (1) k , znk) =N (y(1)n | W (1) xkxn+ W (1) tk tn+ µ (1) k , Ψ (1) k ) znk p(y(2)n | tn, xn; W (2) xk, W (2) tk , µ (2) k , Ψ (2) k , znk) =N (y(2)n | W (2) xkxn+ W (2) tk tn+ µ (2) k , Ψ (2) k ) znk p(znk; πk) = Multinomial(znk| πk)2
ここでN は,以下のような正規分布の確率密度関数である. N (x | µ, Ψ) = √ 1 | 2πΨ |exp [ −1 2{(x − µ) T Ψ−1(x− µ)} ] ゆえに,提案モデルにおいて,観測変数y(1)n およびyn(2)は, 影響を取り除きたい変数xnおよび潜在変数tnによって,平均 が定まる正規分布に従う変数であるといえる.また,πkは混 合係数で0≤ πk≤ 1,p(znk; πk) = πk, ∑K k=1πk= 1である.
3.2
EM
アルゴリズムによる最尤推定
{y, z}に対する対数尤度関数は以下で表される. ln p(y, z| Θ) = N ∑ n=1 K ∑ k=1 znkln πkN (yn| Wxkxn+ µk, Ck) ここでΘはパラメータの組であり,Θ = [πk, µk, ΨkWxk, Wtk] である.また,負担率rnkを以下のように定義する. rnk= πkN (yn| Wxkxn+ µk, Ck) ∑K j=1πjN (yn| Wxjxn+ µj, Cj) EMアルゴリズムによって,パラメータの最尤解を求める. Eステップでは現在のパラメータΘcurを用いて負担率を求め, 完全データに対する対数尤度の期待値Q(Θ, Θcur)を計算する. Q(Θ, Θcur) = N ∑ n=1 K ∑ k=1 rnk{ln πk+N (yn| Wxkxn+ µk, Ck)} Mステップでは求めた対数尤度の期待値を最大化するパラメー タΘnewを求める. Θnew= argmax Θ Q(Θ, Θcur) 以上より,パラメータの更新式は以下のようになる. rnk= πkN (yn| Wxkxn+ µk, Ck) ∑K j=1πjN (yn| Wxjxn+ µj, Cj) µk= ∑N n=1rnk(yn− Wxkxn) ∑N n=1rnk πk= 1 N N ∑ n=1 rnk Wxk= ( N ∑ n=1 rnky˜nkx˜Tnk ) ( N ∑ n=1 rnkx˜nkx˜Tnk )−1 Wtk= Uk(Λk− Ψk) 1 2R Ψk= Sk− WtkWtkT ただし,y˜nk= yn− ∑N n=1rnkyn ∑N n=1rnk , ˜xnk= xn− ∑N n=1rnkxn ∑N n=1rnk , Sk= ∑N n=1rnk( ˜ynk−W∑xkx˜nk)( ˜ynk−Wxkx˜nk)T N n=1rnk であり,UkはSk の固有ベクトルからなる行列,ΛkはSkの固有値を対角成分 に持つ対角行列.Rは任意の直交行列である. また,提案手法において,2.3章と同様に,偏正準相関関係 を求めたい二つの変数をXt,Y (m) t−1 とし,影響を除きたい変 数をXt(m)−1 とすることで,同一のグレンジャー因果で説明可 能なデータ集合を取り出すことができる. 図3: 混合確率的偏正準相関分析のグラフィカルモデル4.
人工データを用いた評価実験
本節では,提案モデルとパラメータ推定のためのEMアル ゴリズムが期待通りに求まるかを検討するため,提案モデルに 適当なパラメータを代入して人工データを生成し,その人工 データからクラスタを推定することを試みる.人工データを生 成するために用いる混合クラスタ数K = 3とする.また,観 測データをy(1)∈ R2,y(2)∈ R3,x∈ R3とする.今回の実 験では,正則化項ηC,ηWxはともに0とした.クラスタ毎の パラメータµ(1),µ(2),W(1) x ,Wx(2),Wt(1),W (2) t ,Ψ (1),Ψ(2),混 合係数πを表1の値に設定し,生成モデルからデータを作成 する. 生成したデータを提案手法である混合確率的偏正準相関分析 および朴らの手法[5]によって解析する.生成モデルから1000 点のデータを作成し,それぞれの手法でクラスタリングを行う までを1試行とする.ここで,あるデータ点xnに対する真の クラスタをk∗n,解析によって得られたクラスタをknとする. 提案手法と既存手法によって得られた結果を比較するために, 各試行における真のクラスタを再現することができなかった割 合,クラスタ非再現率を求める.クラスタ非再現率を定義する ために,以下のように多数派クラスタk(k)ˆ を定義する. ˆ k(k) = argmax j |C k∩ Cj∗| ここで,Ck∗={xn| k∗n= k}は,ある真のクラスタに所属す るデータ集合,Ck{xn| kn= k}は,解析によって得られたあ るクラスタに所属するデータ集合である.多数派クラスタは, 解析によって得られたクラスタに,最も多く存在していた真の クラスタを表している.多数派クラスタを用いて,クラスタ非 再現率を以下のように求める. 1. あるデータ点xnに対して,多数派クラスタˆknと真のク ラスタk∗nが一致していないとき,このデータ点を誤り点 とする. 2. 全てのデータ点に対して誤り点を求め,全てのデータ点 に対して誤り点が占める割合をクラスタ非再現率とする. クラスタ非再現率は,提案手法と真のクラスタの対応関係を多 数派クラスタという形で求めた上での,多数派クラスタに当て はまらない点の割合といえる. それぞれの手法で1000試行の解析を行ったときのクラスタ 非再現率のヒストグラムを図4に示す.図より,朴らの手法 [5]と比較して,提案手法によって正しいクラスタを推定でき ていることが確認できる.3
表1: 人工データのパラメータ(生成モデルからサンプリングした人工データを用いた実験) k µ(1) µ(2) W(1) x Wx(2) 1 ( 0.1 0.2 ) 0.20.3 0.4 (0.1 0.2 0.3 0.4 0.2 0.3 0.4 0.5 ) 0.10.2 0.20.3 0.30.4 0.40.5 0.3 0.4 0.5 0.6 2 ( −0.5 3.5 ) 1.24.3 0.4 (0.5 0.4 −0.3 0.4 −0.1 0.3 0.3 −0.5 ) 0.15.2 −0.30.2 −0.30.4 4.51.4 −0.3 0.4 0.5 −0.6 3 ( 0.1 4.3 ) −4.32.2 2.4 (−1.1 0.7 −0.8 0.4 0.2 2.3 0.4 0.5 ) −0.4 0.7−0.2 0.3 −0.40.3 −0.40.5 0.3 0.1 0.1 0.9 k Wt(1) Wt(2) Ψ(1) Ψ(2) π 1 ( 0.1 0.2 0.2 0.3 ) 0.10.2 0.20.3 0.3 0.4 (0.05 0.08 0.08 0.13 ) 0.210.24 0.240.29 0.320.40 0.32 0.40 0.57 0.33 2 ( 1.1 0.2 0.2 2.3 ) 4.13.2 0.20.3 −2.3 0.4 (0.40 0.18 0.18 0.13 ) 0.890.44 0.440.29 0.800.44 0.80 0.44 0.77 0.33 3 ( −1.0 0.2 0.2 −0.5 ) 0.00.9 −0.6−0.9 0.7 0.1 (1.0 0.0 0.0 0.09 ) 2.931.84 1.841.69 2.931.96 2.93 1.96 6.9 0.34 図4: クラスタ非再現率のヒストグラム
5.
おわりに
本研究では,複数の偏正準相関関係が混在したデータに対 して,それぞれの偏正準相関関係ごとにデータを分割する手法 を提案した.また,提案手法を人工データに適用し,既存手法 [5]と比較して,より正しくクラスタを推定できることを確か めた.紙幅の関係で紹介しなかったが,2者間のコミュニケー ションに関する実データ解析も行っている.本研究の発展とし て,クラスタ数を自動で推定するモデルなどが考えられる.参考文献
[1] B. Raja Rao. Partial canonical correlations. Traba-jos de estadistica y de investigaci´on operativa, Vol. 20, No. 2, pp. 211–219, 1969.
[2] Andr´e Fujita, Joao Ricardo Sato, Kaname Kojima, Lu-ciana Rodrigues Gomes, Masao Nagasaki, Mari Cleide Sogayar, and Satoru Miyano. Identification of granger causality between gene sets. Journal of Bioinformatics
and Computational Biology, Vol. 8, No. 04, pp. 679– 701, 2010.
[3] Yusuke Mukuta and Tatsuya Harada. Probabilistic Partial Canonical Correlation Analysis. In Proceed-ings of the 31st International Conference on Machine Learning, pp. 1449–1457, 2014. [4] ビショップCM. パターン認識と機械学習. シュプリン ガージャパン, 2008. [5] 朴鍾範,森裕紀,浅田稔. Granger causalityに基づく多次 元時系列データ からのcausal pattern抽出. システム制 御情報学会研究発表講演会講演論文集, 2014. Vol.DVD-ROM.
[6] Clive W.J. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: Journal of the Econometric Society, pp. 424–438, 1969.
[7] Harold Hotelling. Relations between two sets of vari-ates. Biometrika, pp. 321–377, 1936.
[8] David Weenink. Canonical correlation analysis. In Proceedings of the Institute of Phonetic Sciences of the University of Amsterdam, Vol. 25, pp. 81–99, 2003. [9] Yuya Yamashita, Tatsuya Harada, and Yasuo
Ku-niyoshi. Causal flow. In IEEE International Confer-ence on Multimedia and Expo, pp. 1–6. IEEE, 2011. [10] Jo˜ao R. Sato, Andr´e Fujita, Elisson F. Cardoso,
Car-los E. Thomaz, Michael J. Brammer, and Edson Amaro Jr. Analyzing the connectivity between re-gions of interest: an approach based on cluster granger causality for fmri data analysis. Neuroimage, Vol. 52, No. 4, pp. 1444–1455, 2010.