正準相関分析を使うと、2つの多次元データ同士の関連性を分析できるらしい。
面白そうなので試してみた。ちなみに正準相関はsklearn.cross_decomposition.CCAで使える。正準相関自体の解説はほとんどしないので、文中のリンクを参考にして欲しい*1。
目次
スポンサーリンク
一応概要だけ
代表的な多変量解析の手法（といって良いのかどうか少し悩むけど）として、主成分分析や重回帰分析が存在する。
- 主成分分析：一つの多変量データを直交するより少ない変数に縮約する
- 重回帰分析：一つの多変量データを一つの単変量データに変換する
主成分分析にしろ重回帰分析にしろ、変換の係数だったり行列だったりを求めてそれで変換するのが実際にやることである。
さて、正準相関は上の流れで説明すると、
- 正準相関分析：二つの多変量データをそれぞれ直交するより少ない変数に縮約して、かつ二つの変換されたデータの間で相関を最大化する
という目的の分析である。主成分分析と重回帰を混ぜた感じ。
気づいた人もいると思うけど、多変量vs多変量のデータでどちらかを単変量に分解して個別に重回帰で解くことも可能である。それに対するメリットとしては、
- 個別に重回帰するより全体の構造みたいなものを捉えられる可能性がある
- 個別に重回帰すると係数の数が全体でとても多くなるので解釈が面倒くさいが、一度次元を下げて直交した空間に持っていくことでそこが楽になる
というあたりがあり、要するに解釈性がいいということ。
この説明でもよくわからん、という人は、ニューラルネットのオートエンコーダーとか思い浮かべていただくとかえってわかりやすいかもしれない。
ノイズに埋もれた波形を取り出す
参考URLの通りにやることにする。
単一の信号源に複数のプローブを当てていて、それぞれに独立のノイズが乗って信号が埋もれてしまった・・・みたいな状況から元の信号を取り出そうとしているらしい。脳波計測とかで使えるのだろうか？
参考URL：https://www.jstage.jst.go.jp/article/jnns/20/2/20_62/_pdf
とりあえずこのようなコードを書き、
# coding: UTF-8 import numpy as np import matplotlib.pyplot as plt from sklearn.cross_decomposition import CCA def plot_wave(data, filename): fig, [ax1,ax2] = plt.subplots(2,1,figsize=(16,9)) ax1.plot(data[:,0], color="b") ax2.plot(data[:,1], color="r") plt.savefig(filename) def gen_pulse_data(): common_pulse = np.array([-1]*50 + [0]*50 + [1]*50 + [0]*50 + [-1]*50 + [0]*50 + [1]*50 + [0]*50, dtype=np.float64) common_pulse += (np.random.random(common_pulse.shape) - 0.5)*0.1 noise1 = (np.random.random(common_pulse.shape) - 0.5)*50 noise2 = (np.random.random(common_pulse.shape) - 0.5)*50 data1 = np.vstack([common_pulse + noise1, common_pulse - noise1]).T data2 = np.vstack([common_pulse + noise2, common_pulse - noise2]).T return data1, data2 def main(): X1, X2 = gen_pulse_data() plot_wave(X1, "X1.png") plot_wave(X2, "X2.png") cca = CCA(n_components=2) cca.fit(X1, X2) Y1 = cca.transform(X1) Y2 = cca.transform(X2) plot_wave(Y1, "Y1.png") plot_wave(Y2, "Y2.png") if __name__ == "__main__": main()
実行する。
まずは元の信号。
X1
X2
わかる訳ねえな、という感じ。
CCAでX1とX2の相関が最大になるような変換を計算し、その変換に基いてX1を変換したものをY1とすると、
Y1
こんな感じになった。Y2も同じようなものなので省略。ここでは2次元出しているが、元のパルス信号が1次元なので2次元目（下）はノイズだけ出ている。
まあ、上手く動いているのではないだろうか。
もうちょっとデータ分析っぽいことをしてみる
如何せん、「ノイズに埋もれた信号を取り出せる！」というだけでは、データ分析っぽくなくて（個人的には）面白くない。正準相関自体はもっと色々なことができるはず。
ここで足を引っ張るのは「正準相関向きのサンプルデータが見つからない」ということである。
正準相関向きのデータを探すのは困難
2つの多次元データが対応しているようなデータで、適当にわかりやすいものがあれば良いのだが・・・なかなか良いデータがない。上に挙げた解説論文でも、「知名度は低い」とか書かれちゃってるし、正準相関自体、ニッチな感じがする。そこが素敵なのだが。
一応、ネット上にある解説例だと、
医学の分野で、肝機能の検査値（複数）と腎機能の検査値（複数）の対応を見るとか、
中学生の体格（身長、体重、座高とか）と運動能力（50m走、走り幅跳びとか）の対応を見るとか、
そういう感じのことをやっているのだが、この手のデータを探してくるのがまず面倒くさいし、見つけてもプログラムに流し込めるようにするまでがまた苦行だろうな、ということは容易に想像できるのである。
この点で悩んで、この記事も一週間くらい出すか出さないか迷ってたんだけど、やることにした。ただ、結局良いデータは見つからなかったので、それっぽくでっちあげることにした。
作成したデータ
ある架空の中学校で集計したという設定の、20人の生徒のデータである。「学外での勉強や取り組み」と「学校の成績」が対応付いている。
「学外での勉強や取り組み」には、
- 一ヶ月に何冊読書するか
- 一年に何回博物館に行くか
- 毎週何日塾に通っているか
- 毎日何時間自習しているか
の4つの変数がある。一方、「学校の成績」は、
- 国語
- 数学
- 社会
- 理科
- 英語
の5つの科目があり、5段階評価で成績が付く。
本来であれば適当に線形モデルでも作ってあげて数字を作るべきところだが、面倒くさいので私の想像で適当に埋めた（ツッコミポイント）。
一応、次のような方針を考え、それに沿ったデータになるようにでっちあげた・・・つもり。
- 読書量と国語の成績は比例する
- 博物館に行った回数と社会、理科の成績は比例する
- 塾に通う頻度、自習時間は成績全体に影響を及ぼす
よって、こういう結果が出てくるか、という勝負になる。
実験と結果
こういうプログラムを書いた。
# coding: UTF-8 import numpy as np from scipy.stats import pearsonr from sklearn.cross_decomposition import CCA def gen_data(): # X1: # 毎月何冊の本を読むか, # 一年に何回博物館に行くか, # 塾に週何日通うか, # 毎日何時間自習するか X1 = np.array([[1,0,2,1], [3,2,4,2], [0,0,2,0], [9,4,2,1], [1,1,3,1], [8,1,6,3], [0,9,7,8], [2,2,4,1], [5,0,0,1], [2,0,4,0], [0,0,7,8], [4,4,2,2], [5,1,2,1], [1,1,5,2], [8,6,2,1], [0,0,0,1], [6,1,3,1], [2,0,3,1], [4,8,5,3], [5,0,1,1]]) # X2: # 国語,数学,社会,理科,英語の成績 X2 = np.array([[3,3,3,3,3], [4,3,4,4,5], [2,2,3,3,2], [5,4,3,3,3], [3,3,4,4,4], [5,5,5,4,5], [3,5,5,4,5], [4,4,4,5,3], [5,3,3,3,3], [3,4,3,4,3], [5,5,4,5,5], [4,4,5,5,3], [4,3,3,3,3], [4,4,5,4,5], [5,3,5,5,3], [2,2,2,1,2], [5,3,4,4,4], [3,4,3,4,3], [5,5,5,5,5], [5,3,3,3,3]]) return X1, X2 def main(): X1, X2 = gen_data() cca = CCA(n_components=4) cca.fit(X1, X2) print("Correlation Coefficient") for i in range(4): print("{0}:{1:.4f}".format(i, pearsonr(cca.x_scores_[:,i], cca.y_scores_[:,i])[0])) print("") np.set_printoptions(formatter={'float': '{: 0.4f}'.format}) print("X1 loadings") print(cca.x_loadings_.T) print("") print("X2 loadings") print(cca.y_loadings_.T) if __name__ == "__main__": main()
「学外での勉強や取り組み」=X1と「学校の成績」=X2を4次元の空間上に写像して相関を最大化する、という問題を解かせる。軸同士は直交していて無相関なので、写像したデータの各軸の値同士の相関だけ見てやれば良い。写像したデータは、cca.x(or y)_scores_かcca.transform(X1(or X2))で取得できる*2。
あとはX1とX2の各成分が、写像先の各軸にどれだけ寄与しているかがわかれば良い。そのためにはcca.x(or y)_loadings_を見る。転置した方が見やすいのでそうしている。
こうして見ると、PCAに似ている。実際、CCAとPCAは親戚らしい。ま、あまり理論的な話に深入りしてもボロが出るので、これくらいにしておく。
さて、結果はこのようになった。
Correlation Coefficient 0:0.9558 1:0.8978 2:0.5980 3:0.2927 X1 loadings [[-0.4224 0.4244 1.0047 0.7353] [ 0.9326 0.4315 0.2511 0.3450] [ 0.2350 0.9269 -0.1714 -0.2514] [-0.4357 0.4062 -0.0634 0.8007]] X2 loadings [[-0.0558 0.7802 0.6574 0.6206 0.8044] [ 0.9392 0.5254 0.4702 0.3366 0.4271] [-0.2184 -0.0705 0.9353 0.5442 -0.3528] [-0.4310 0.0676 -0.2136 -0.8751 0.0008]]
まず見るべきはCorrelation Coefficientで、写像先の空間の軸にどれだけ相関（＝やった意味）があるかを示している。0,1,2次元目はまあまあ強い相関だが、3次元目は相関係数0.3じゃ大した意味はなさそうだな、という風に解釈しておく。
次にX1 loadingsとX2 loadingsを見る。X1 loadingsは4*4、X2 loadingsは4*5で、つまり行が写像先の軸、列が元の空間の軸に対応するように表示している。
X1 loadingsの各行を見ていくと、
- 1行目
塾と自習に熱心
- 2行目
読書
- 3行目
博物館
- 4行目
自習と博物館だけ？
なんとか解釈できる。数字がでかいところだけ重視するのがこつ。X2 loadingsも同様にやると、
- 1行目
国語以外のすべて。国語にはほぼ中立。特に強いのは英語
- 2行目
国語。他もそれなりに
- 3行目
社会と理科
- 4行目
理科にとてもネガティブ。全体的にネガティブな感じ
ここまで出揃えば後はなんとかなる。このデータを作った方針を再掲する。
- 読書量と国語の成績は比例する
- 博物館に行った回数と社会、理科の成績は比例する
- 塾に通う頻度、自習時間は成績全体に影響を及ぼす
0次元目は「塾に通う頻度、自習時間は成績全体に影響を及ぼす」とに、1次元目は「読書量と国語の成績は比例する」に、2次元目は「博物館に行った回数と社会、理科の成績は比例する」に対応していることがわかり、まあ妥当な結果なんじゃないの、という気はする。相関係数の低い3次元目はそこまで重視する必要はない。
今回は先に方針を決めてデータをでっち上げたのであまり感動がないような気もするが、実際はデータにどんな構造があるのかは分析してみないとわからない。その構造を理解する上で正準相関が役に立つことは、上の例でなんとなく理解できた。