[[FrontPage]]
#contents
2011/06/29からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://sage.math.canterbury.ac.nz/home/pub/108/
また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* 第9章-混合ガウス分布のEMアルゴリズムをSageで試す [#mbf01df5]
混合ガウス分布の最尤推定問題をEMアルゴリズムで解いてみます。
ここで紹介するのは、PRMLの図9.8、
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8a.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8b.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8c.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8d.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8e.jpg,center,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.8f.jpg,center,30%);
です。反復を繰り返して収束する様子がよく分かります。
** テストデータ [#b83b4ad9]
テストデータには、Old Faithful間欠泉データをX-Y平面に-2.5から2.5の範囲に正規化したものを使用します。
sageへの入力:
#pre{{
# 混合ガウス分布のEMアルゴリズム
from numpy import loadtxt
data = loadtxt(DATA+"faithful.txt")
# 1カラムはx, 3カラムにはyがセットされている
# これを-2.5, 2.5に正規化する
x = data[:,0]- 3.5
y = (data[:,1]-70)*5/60
# プロット
data_plt = list_plot(zip(x,y))
data_plt.show(aspect_ratio=1, figsize=(3), xmin=-2.5, xmax=2.5, ymin=-2.5, ymax=2.5)
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# ベクトルのガウス分布を定義
def _gauss(v, mu, sigma):
d = len(v);
sigma_inv = sigma.inverse();
sigma_abs_sqrt = sigma.det().sqrt();
# sage 4.6.2ではtranspose()の代わりにcolumn()を使用する
val = -(v - mu) * sigma_inv * (v - mu).transpose()/2;
a = (2*pi)**(-d/2) * sigma_abs_sqrt**-0.5
return a * e**val[0];
}}
sageへの入力:
#pre{{
# 対数尤度
def ln_p():
sum = 0.0
for n in range(N):
sum += ln(pi_N[n].sum())
return sum
}}
sageへの入力:
#pre{{
# (x - u)(x - u)^Tを計算
def covMatrix(x, u):
d = x - u
return matrix(RDF, [[d[i]*d[j] for j in range(D)] for i in range(D)])
}}
sageへの入力:
#pre{{
# 定数の初期化
D = 2
K = 2
N = len(x)
beta = 0.5
X = matrix(RDF, zip(x, y))
ident = identity_matrix(2)
# アクションのインデックス
action_idx = {0:0, 1:1, 5:5, 20:20}
}}
sageへの入力:
#pre{{
# アクション関数
def action_f(mu_k, sigma_k, gamma):
pt_plt = Graphics()
for n in range(N):
pt_plt += point(X[n], rgbcolor=(gamma[n][1], 0, gamma[n][0]))
r_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu_k[1], sigma_k[1]), [x, -2.5, 2.5], [y, -2.5, 2.5],
contours = 1, cmap=['red'], fill=False)
b_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu_k[0], sigma_k[0]), [x, -2.5, 2.5], [y, -2.5, 2.5],
contours = 1, cmap=['blue'], fill=False)
(r_cnt_plt + b_cnt_plt + pt_plt).show(aspect_ratio=1, figsize=(3), xmin=-2.5, xmax=2.5, ymin=-2.5, ymax=2.5)
}}
** EMアルゴリズム [#o9672258]
混合ガウス分布のEMアルゴリズムは、以下の通りです。
+ 平均\(\mu_k\)、 分散\(\Sigma_k\)、混合係数\(\pi_k\)を初期化する
+ Eステップ：式(9.23)を使って負担率を計算する
$$
\gamma(z_{nk}) = \frac{\pi_k \mathcal{N}(x_n|\mu_k, \Sigma_k)}{\sum_{j=1}^N \pi_j \mathcal{N}(x_n|\mu_j, \Sigma_j)}
$$
+ Mステップ：上記の負担率を使って式(9.24), (9.25), (9.26), (9.27)のパラメータ値を再計算する
$$
\mu_k^{new} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma(z_{nk})x_n
$$
$$
\Sigma_k^{new} = \frac{1}{N_k}\sum_{n=1}^{N} \gamma(z_{nk}) (x_n - \mu_k^{new}) (x_n - \mu_k^{new})^T
$$
$$
\pi_k^{new} = \frac{N_k}{N}
$$
$$
N_k = \sum_{n=1}^{N} \gamma(z_{nk})
$$
+
対数尤度：式(9.28)を計算し、対数尤度が収束するまでステップ２に戻る
$$
\ln p(X|\mu, \Sigma, \pi) = \sum_{n=1}^{N} ln \left \{ \sum_{k=1}^{K} \pi_k \mathcal{N}(x_n|\mu_k, \Sigma_k) \right \}
$$
sageへの入力:
#pre{{
# EMアルゴリズム
def _EM(pi_k, mu_k, sigma_k):
# 対数尤度の初期値
o_lnP = ln_p()
diff = 1
for l in range(21):
# Eステップ
pi_N = matrix(RDF, [[pi_k[k]*_gauss(X[n], mu_k[k], sigma_k[k]) for k in range(K)] for n in range(N)])
gamma = matrix(RDF, N, K)
for n in range(N):
gamma.set_row(n, pi_N[n]/(pi_N[n].sum()))
# 最初のEステップの状態
if l == 0:
action_f(mu_k, sigma_k, gamma)
# Mステップ
N_k = [gamma.column(k).sum() for k in range(K)]
for k in range(K):
mu_k[k] = 0
for n in range(N):
mu_k[k] += gamma[n][k]*X[n]
mu_k[k] /= N_k[k]
sigma_k[k] = 0
for n in range(N):
sigma_k[k] = sigma_k[k] + gamma[n][k]*covMatrix(X[n], mu_k[k])
sigma_k[k] /= N_k[k]
pi_k[k] = N_k[k]/N
# 新しい対数尤度
lnP = ln_p()
diff = abs(lnP - o_lnP)
o_lnP = lnP
# 図化
if action_idx.has_key(l):
action_f(mu_k, sigma_k, gamma)
}}
** 計算結果 [#h2dca0ef]
平均の初期値を\(\mu_1 = (-1.5, 0.5), \mu_2 = (1.5, -0.5)\)、 分散の初期値を0.5として、
計算した結果を以下に示します。
図9.8に近い結果を得ることができました。
sageへの入力:
#pre{{
# mu=(-1.5, 0.5), (1.5, -0.5)で計算
pi_k = [0.5, 0.5]
mu_k = [vector([-1.5, 0.5]), vector([1.5, -0.5])]
sigma_k = [matrix(RDF, D, D, beta*ident), matrix(RDF, D, D, beta*ident)]
_EM(pi_k, mu_k, sigma_k)
}}
&ref(sage0-1.png);
&ref(sage1.png);
&ref(sage2.png);
&ref(sage3.png);
&ref(sage4.png);
** μの初期値に依存 [#n91c2979]
平均の初期値が図9.8と少しずれていたので、\(\mu_1 = (-1.5, 1.0), \mu_2 = (1.5, -1.0)\)
としたところ、途中まで上手く計算していたのですが、L=20回目がまったく異なる結果になってしまい
ました。
これはもしかすると、図9.7
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure9.7.jpg,center,40%);
で説明されている特異性の例かも知れない。
sageへの入力:
#pre{{
# μの初期値に結果が大きく依存する
# mu=(-1.5, 1.0), (1.5, -1.0)で計算
pi_k = [0.5, 0.5]
mu_k = [vector([-1.5, 1.0]), vector([1.5, -1.0])]
sigma_k = [matrix(RDF, D, D, beta*ident), matrix(RDF, D, D, beta*ident)]
_EM(pi_k, mu_k, sigma_k)
}}
&ref(sage0-2.png);
&ref(sage1-1.png);
&ref(sage2-1.png);
&ref(sage3-1.png);
&ref(sage4-1.png);
** コメント [#f724fe13]
#vote(おもしろかった[16],そうでもない[1],わかりずらい[0])
#vote(おもしろかった[17],そうでもない[1],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
