sage/PRML - 混合ガウス分布のEMアルゴリズム
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/29からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/108/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第9章-混合ガウス分布のEMアルゴリズムをSageで試す [#mbf0...
混合ガウス分布の最尤推定問題をEMアルゴリズムで解いてみま...
ここで紹介するのは、PRMLの図9.8、
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
です。反復を繰り返して収束する様子がよく分かります。
** テストデータ [#b83b4ad9]
テストデータには、Old Faithful間欠泉データをX-Y平面に-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, xma...
}}
&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...
}}
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, g...
r_cnt_plt = contour_plot(lambda x, y : _gauss(vector(...
contours = 1, cmap=['red'], fill=False)
b_cnt_plt = contour_plot(lambda x, y : _gauss(vector(...
contours = 1, cmap=['blue'], fill=False)
(r_cnt_plt + b_cnt_plt + pt_plt).show(aspect_ratio=1,...
}}
** EMアルゴリズム [#o9672258]
混合ガウス分布のEMアルゴリズムは、以下の通りです。
+ 平均\(\mu_k\)、 分散\(\Sigma_k\)、混合係数\(\pi_k\)を初...
+ Eステップ：式(9.23)を使って負担率を計算する
$$
\gamma(z_{nk}) = \frac{\pi_k \mathcal{N}(x_n|\mu_k, \Sigm...
$$
+ Mステップ：上記の負担率を使って式(9.24), (9.25), (9.26)...
$$
\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...
$$
$$
\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 \{ \s...
$$
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]...
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]*co...
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)\...
計算した結果を以下に示します。
図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, ...
_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,...
としたところ、途中まで上手く計算していたのですが、L=20回...
ました。
これはもしかすると、図9.7
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
で説明されている特異性の例かも知れない。
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, ...
_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(おもしろかった[17],そうでもない[1],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/29からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://sage.math.canterbury.ac.nz/home/pub/108/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* 第9章-混合ガウス分布のEMアルゴリズムをSageで試す [#mbf0...
混合ガウス分布の最尤推定問題をEMアルゴリズムで解いてみま...
ここで紹介するのは、PRMLの図9.8、
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
です。反復を繰り返して収束する様子がよく分かります。
** テストデータ [#b83b4ad9]
テストデータには、Old Faithful間欠泉データをX-Y平面に-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, xma...
}}
&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...
}}
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, g...
r_cnt_plt = contour_plot(lambda x, y : _gauss(vector(...
contours = 1, cmap=['red'], fill=False)
b_cnt_plt = contour_plot(lambda x, y : _gauss(vector(...
contours = 1, cmap=['blue'], fill=False)
(r_cnt_plt + b_cnt_plt + pt_plt).show(aspect_ratio=1,...
}}
** EMアルゴリズム [#o9672258]
混合ガウス分布のEMアルゴリズムは、以下の通りです。
+ 平均\(\mu_k\)、 分散\(\Sigma_k\)、混合係数\(\pi_k\)を初...
+ Eステップ：式(9.23)を使って負担率を計算する
$$
\gamma(z_{nk}) = \frac{\pi_k \mathcal{N}(x_n|\mu_k, \Sigm...
$$
+ Mステップ：上記の負担率を使って式(9.24), (9.25), (9.26)...
$$
\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...
$$
$$
\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 \{ \s...
$$
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]...
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]*co...
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)\...
計算した結果を以下に示します。
図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, ...
_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,...
としたところ、途中まで上手く計算していたのですが、L=20回...
ました。
これはもしかすると、図9.7
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
で説明されている特異性の例かも知れない。
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, ...
_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(おもしろかった[17],そうでもない[1],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
