狠狠撸

Numpy/Scipyで
独立成分分析
2011年8月28日
第1回Tokyo.Scipy
@sfchaos

目次

1. 自己紹介
2. 独立成分分析
3. Numpy/Scipyを用いた実装
4. まとめ

1

目次

1. 自己紹介
4. まとめ

2

1.1 プロフィール
? TwitterID: @sfchaos

? 職業：コンサルタント
? 金融工学のモデル構築?データ解析
? 最近，大規模データ解析の企画に着手
(Hadoop, Mahout, CEP等)

3

1.2 RとPythonと私

? R歴：4年
? 「Rパッケージガイドブック」(東京図書，2011年4月)執筆．
? bigmemoryパッケージ(大規模データの管理?分析)
? RTisean/tseriesChaosパッケージ(非線形(カオス)時系列解析)

? Tokyo.R 翻訳プロジェクト Wiki
? Tokyo.Rの方々と各種マニュアルの翻訳

4

? Python歴：薄く1年(Numpy/Scipy歴1ヶ月)

? Rでは大規模データの扱いが大変なので，代替となる手
段がないか模索中
→ PythonでのHPCに興味があります．

? 今回は，Numpy/Scipyの勉強のために発表させていた
だきます(独立成分分析自体も初心者)．

5

目次

1. 自己紹介
4. まとめ

6

2.1 独立成分分析の概要
? 未知音源分離問題(Blind Source Separation, BSS).
　観測信号から元の信号の成分を取り出す問題．カクテ
ルパーティー効果など．
? 脳磁波の解析，音声処理などで使用されている(らし
い)．
独立成分分析では，
観測信号のみから原信号を推定

原信号1 観測信号1

原信号i
重ね合わせ観測信号i

原信号N 観測信号N
7

2.2 無相関性と独立性

? 日常用語では無相関も独立も同じ意味で使われること
が多い
? しかし，両者の数学的な定義は異なり，独立?無相関

確率変数 x1と x2 が無相関
無相関
E[ x1 x 2 ] = E[ x1 ]E[ x 2 ]
独立
確率変数 x1と x2 が独立
p1, 2 ( x1 , x2 ) = p1 ( x1 ) p2 ( x2 )

8

? 原信号がそれぞれ正規分布に従う場合は，「無相関＝
独立」となり，主成分分析を行えば元信号を取り出せる．
? 主成分分析は，各主成分が正規分布に従うと仮定して
おり，正規分布に従う変数を線形結合した変数も正規
分布に従うことに注意(PRML 第12章，2章)．

9

? しかし，原信号が正規分布に従わない場合，線形結合
した変数同士が無相関だが独立ではないことがある．
観測信号は無相関だが
原信号は無相関かつ独立
独立ではない

? s1 ? 原信号の重ね合わせ ? x1 ? ? s1 ?
? ?
?s ? ? ?=
?x ? A? ?
?s ?
? 2? ? 2 3? ? 2? ? 2?
A=?
? 2 1?
?
? ? p( x1 , x 2 ) ≠ p1 ( x1 ) p2 ( x 2 )
10

? このような場合，主成分分析を行っても無相関ではある
が独立ではない原信号が得られてしまう．

無相関だが独立ではない

12

? このような観測信号から，独立な原信号の成分を抽出
することが独立成分分析の目的．
観測信号原信号

独立
? x1 ? ? x1 ?
? ?
?x ? W? ?
?x ?
? 2? W ? 2?
13

2.3 独立成分の抽出

? 目標：この2つの方向を取り出すこと
? 目で見れば一目瞭然だが，機械的に抽出したい

観測信号

14

2.3.1 独立性
? 独立な方向を特徴付ける指標を定義したい．
? 実は，分布の非正規性を最大化することが独立成分を
求めることと等価であることが分かっている(きちんとし
た議論は「詳解独立成分分析」を参照)．
? 非正規性を特徴付ける指標はいろいろある．
概要
尖度(の絶対値) ? 4次のモーメント
? 正規分布では3次以上のキュムラントはゼロであり，尖度が大き
いほど正規分布から離れていることになる
ネゲントロピー ? 正規分布は平均と分散が一定の分布の中で最も情報量が大き
い(ランダムである)ことに着目
? 正規分布と対象とする分布の情報量の差をネゲントロピーとして
定義
? ネゲントロピーが大きい分布ほど正規分布から離れていることに
なる
他にも，相互情報量などがある． 15

? 尖度(4次のモーメント)
( x) = E[ x ] ? 3{E[ x ]} ? 4 E[ x ]E[ x] + 12 E[ x ]{E[ x]} ? 6{E[ x]}
4 2 2 3 2 2 4
κ4

? ネゲントロピー
J ( y ) = H ( y gauss ) ? H ( y ),
H ( y ) = ? ∫ p ( y ) log p( y ) dy p( y ) ：分布の確率密度関数

ネゲントロピーを計算するためには確率密度関数を近似する必要が
あるので，後ほど近似関数を導入する．

16

2.3.3 独立成分抽出のアイディア
? 求めたい独立成分の個数分だけ，非正規性の指標が
大きい方向を最適化により求める．

wi (0) w (i j )
???

非正規性の指標

max E[ g i ( yi )], i = 1,… , N
yi = wi x

s.t. wi = 1.
17

? これは，観測信号に対して任意の線形変換を行い，独
立な成分を探すという問題であり，自由度が大きいため
大変．

?

18

? そこで，独立?無相関という関係に着目する．
? 観測信号を無相関にしておけば，独立な成分を探索す
るためには直交変換だけを対象とすればO.K.
? 主成分分析を用いて，観測信号を直交する成分に変換
しておく(白色化，whitening, sphering)

白色化

z = Vx

19

白色化

x z = Vx z

独立成分の
抽出
y = Wz

20

2.4 FastICAアルゴリズム
? 収束が速い(3次収束)．
? 最適化手法としてNewton法を使用．
? 独立性の指標は，ネゲントロピーを使用．

max E[ g i ( y i )], i = 1,… , N λ T
yi = wi x
E[ g i ( y i )] + w i w i = max .
s.t . w i = 1. 2

21

λ T
E[ g i ( yi )] + w i w i = max .
2

上式を微分すると，次の方程式の解を求める問題に帰着する．
f ( x i ) = E[ x i g i ( yi )] + λw i = 0.
Newton法を用いて，重みベクトルを求める．
w i( k +1) = w i( k ) ? f ( z i ) / f ' ( z i )
E ( z i g ( w i( k ) )) + λw i( k )
= w i( k ) ? T
E ( g ' ( w i( k ) z i )) + λ
T
∝ E ( z i g ( w i( k ) ) w i( k ) ) ? E ( g ' ( w i( k ) z i ))

ネゲントロピーは，以下の式で近似することが多い．
g1 ( y ) = tanh(ay ), g 2 ( y ) = y exp(? y 2 / 2), g 3 ( y ) = y 3 ,
22

FastICA
p : 推定する独立成分の個数, m :独立成分のカウンタ(初期値=1)
1.(中心化)　データを中心が原点になるように変換する：
1 n
x i' = x i ? ∑ x j
'
n j =1
2.(白色化)　　 ' を白色化(主成分分析)して，データ点 i を得る：
x z
i
3. 各独立成分に対するの射影方向　　p の初期値をランダムに与える．w

m=1, 2, ???, p に対して以下を行う．
w
4. 次式に従って　　 p を更新する：
w p ← E[ z i g ( w T ) w p ] ? E[ g ' ( w T ) z i ]
p p
p ?1
5. 次式に従ってを正規直交化する： w p ? ∑ (w T w p )w j
j
j =1
wp ← p ?1
w p ? ∑ (w T w p )w j
j
j =1

6. w pが収束していない場合は，step.4に戻る．
7. mを1つインクリメントし，の場合はstep.4に戻る．
23

目次

1. 自己紹介
4. まとめ

24


? Pythonには，MDPやscikitsなど独立成分分析を行う
パッケージが既にある．
? ここでは，Numpy/Scipyを用いてFastICAアルゴリズム
を実装する．

25

3.1 使用するNumpy/Scipyの機能?関数
? numpy.arrayオブジェクト等を使用して，行列演算，乱数
生成等にNumpy, Scipyを使用する．

機能関数
特異値分解 linalg.svd関数(Numpy linglg)
行列演算 (主成分分析の中で使用)
ベクトルの正規化 linalg.norm関数(Numpy )
乱数生成乱数生成 random.normal関数(Numpy)

26

3.2　ソースコード
import numpy as np
import scipy as sp
import pylab
from scipy import linalg
# FastICA
def fastICA(x, n_comp, alpha=1.0, maxit=200, tol=1e-4):
n, p = x.shape
# step.1 観測信号の中心化
x = x - x.mean(axis=0)
# step.2 白色化(主成分分析)
v = np.dot(x.T, x)/n
u, sigma, v = linalg.svd(v) # 特異値分解
D = np.diag(1/np.sqrt(sigma))
　　V = np.dot(D, u.T) # 白色化行列の生成
V = V[0:n_comp, 0:p] # 成分の個数
z = np.dot(V, x.T).T # 白色化を実行してベクトルzを得る
# step.3 射影ベクトルの初期化
w_init = np.random.normal(loc=0, scale=1, size=(n_comp, n_comp))
27

# step.4-7 減次法によるFastICAの実行
w = ica_deflation(z, n_comp, alpha, maxit, tol, w_init)
# 復元信号
s = np.dot(w, z.T).T
pylab.scatter(s[:, 0], s[:, 1])
pylab.xlabel("$y_{1}$", fontsize=30)
pylab.ylabel("$y_{2}$", fontsize=30)
pylab.savefig("../output/scatter_estimated_sources.png")
pylab.show()
pylab.close()
return s

28

# 減次法によるFastICA
def ica_deflation(x, n_comp, alpha, maxit, tol, w_init):
n, p = x.shape
# 求める射影ベクトルの結果を格納するオブジェクト
w_res = np.zeros((n_comp, n_comp))
# step.4-7 各独立成分を求める
for i in np.arange(n_comp):
w = w_init[i, ].T
t = np.zeros((n_comp, np.size(w)))
# 既に求められている独立成分と直交するように変換
if i > 0:
k = np.dot(w_res[0:i, ], w)
t = (w_res[0:i, ].T * k).T
w -= t.sum(axis=0)
w /= np.linalg.norm(w)

29

# 誤差，反復回数カウンタの初期化
lim = np.repeat(1000.0, maxit); it = 0
# 独立成分の探索
while lim[it] > tol and it < maxit:
# step.4 wpの更新
wx = np.dot(x, w)
gwx = np.tanh(alpha * wx)
gwx = np.repeat(gwx, n_comp).reshape(np.size(gwx), n_comp)
xgwx = x * gwx
v1 = xgwx.mean(axis=0)
g_wx = alpha * (1 - (np.tanh(alpha * wx))**2)
v2 = np.dot(np.mean(g_wx), w)
w1 = v1 - v2
it += 1
t = np.zeros((n_comp, np.size(w)))
# step.5 Gram-Schmidtの正規直交化
if i > 0:
k = np.dot(w_res[0:i, ], w1)
t = (w_res[0:i, ].T * k).T
w1 -= t.sum(axis=0); w1 /= np.linalg.norm(w1)
# step.6 収束判定
lim[it] = abs(abs(sum(w1 * w)) - 1.0)
w = w1
w_res[i, ] = w
return w_res
30

3.3 独立成分の方向への収束の確認

w1 w2
(4回で収束) (正規直交化するため，
1回で収束する)

①

②
③
④

31

3.4 原信号に関して復元できるとは限らない性質

? 原信号の順番

1番目の
2番目の 1番目の独立成分
独立成分独立成分 2番目の
独立成分

32

? 原信号の符号と強度
? 正負が逆転する可能性がある．
? 射影ベクトルを単位ベクトルとしているため，波形は原信号と
相似形のものが得られ，強度は復元されない可能性がある．

2番目の 1番目の 2番目の 1番目の
独立成分独立成分独立成分独立成分

2番目の 2番目の 1番目の
1番目の
独立成分独立成分独立成分
独立成分
33

目次

1. 自己紹介
4. まとめ

34

4. まとめ
? 独立成分分析は，観測された信号から独立な原信号
の成分を抽出するための多変量解析手法．

? 代表的なアルゴリズムであるFastICAをNumpy/Scipyを
使って実装した．

? Rに比べてNumpy/Scipyが勝っている点は???
→ これを見つけることが今後の課題

35

参考文献

? 「詳解　独立成分分析」

? 朱鷺の杜　Numpy
? NumPy for R (and S-Plus) users
RユーザがPythonの該当関数を調べるのに便利．
http://mathesaurus.sourceforge.net/r-numpy.html

36

ご清聴ありがとうございました

37

狠狠撸

Numpy scipyで独立成分分析

More Related Content

Numpy scipyで独立成分分析