pooneilの脳科学論文コメント

■ EMアルゴリズムの勉強メモ

自由エネルギー原理を理解するためには機械学習での「変分ベイズ」を理解する必要があって、さらにその手前の段階に「EMアルゴリズム」がある。EMアルゴリズムにおいてもKL divergenceを最小化して下界Lを最大化する過程が出てくる(PRMLの9.4章の図9.11-14)。

この図と式を字面を追っていくことはできるけど、シンプルなモデルでじっさいにグラフを書いて理解できるようにしたい。いちばん簡単な例はなんだろうか? Nature BiotechnologyのPrimerの記事で混合二項分布を使ったいい感じにわかりやすいものを見つけた：Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm? Nature Biotechnology, 26(8), 897–899.

ここにある事例を使ってmatlabでグラフ書きながら理解してみることにしよう。

---

潜在変数がない場合 (=>最尤推定)

コインの裏表のデータがある。データXはコインの表の数(列1)と裏の数(列2)を表していて、10回での裏表の数を計算して、それを5回繰り返しす。つまり、データ $X = \{x_1, x_2, ..., x_5\}$ で、たとえば $x_1 = [5\ 5]$ だと、10回中5回表が出たということ。

X = [...
    5 5;
    9 1;
    8 2;
    4 6;
    7 3];

じつはコインには二種類、コインAとコインBがあって、それぞれで表が出る確率が違う。この5回でどっちを使ったかのデータがわかっていて、それをZで示している。つまり、データ $Z = \{z_1, z_2, ..., z_5\}$ で、$z_1 = \{0,1\}$ (A or BでAを使った確率)

Z = [...
    0
    1
    1
    0
    1];

それではcoin Aおよびcoin Bそれぞれでの表の出る確率 $\theta_A, \theta_B$はいくつか？

これは最尤推定で解ける。しかも二項分布の場合は最尤推定を使わなくても、coin Aでは thetaA = 24/30=0.80 , coin Bでは thetaB = 9/20=0.45 という計算だけで済む。

---

潜在変数がある場合 (=>EMアルゴリズム)

ではもし、いつcoin A, coin Bが出たか(Z)が不明の場合にもcoin Aとcoin Bの表が出る確率を推定することはできるだろうか？ つまり入手できるデータはXだけで、Zのほうは入手できない潜在変数という場合でも可能か？というのが問い。

そうすると問題はこのときの尤度

$$p(X,Z) = p(x_1|z_1)p(z_1) * p(x_2|z_2)p(z_2) * ... * p(x_5|z_5)p(z_5)$$

を最大化するような $\theta_A, \theta_B, z_i$ を求めたい、ということになる。

そこで尤度 $p(X,Z)$ を最大化する代わりに、

$p(X,Z) = p(Z|X) * p(X)$ における $p(X)$ を最大化する。

$$\begin{eqnarray} p(x_1) &=& p(x_1, z_1=A) + p(x_1, z_1=B)\\ &=& p(x_1|z_1=A)p(z_1=A) + p(x_1|z_1=B)p(z_1=B)\\ &=& \sum_z{p(x_1|z)p(z)}\\ p(X) &=& p(x_1) * p(x_2) * ... * p(x_5)\\ \end{eqnarray}$$

これのlogを取ったものの最大化を考える。

$$log(p(X)) = \sum_{i}{ log( \sum_{z}{p(xi|z)p(z)} ) }$$

(1) 推定する変数 $\theta$ および $z_i$ の初期化

$\theta_A$, $\theta_B$ : コインA, Bそれぞれの表が出る確率

theta_Aest = 0.60;
theta_Best = 0.50;

$p(z_i=A)$ および $p(z_i=B)$

pziAest = [...
    0.5;
    0.5;
    0.5;
    0.5;
    0.5];
pziBest = 1 - pziAest;

(2) E-step: $p(z_i|x_i)$ を計算する

ベイズの公式より

$$\begin{eqnarray} p(z_i|x_i) &=& p(x_i|z_i)p(z_i) / p(x_i)\\ p(x_i) &=& p(z_i=A,x_i) + p(z_i=B,x_i) \tag{1}\\ p(z_i=A,x_i) &=& p(x_i|z_i=A)p(z_i=A) \tag{2a}\\ p(z_i=B,x_i) &=& p(x_i|z_i=B)p(z_i=B) \tag{2b} \end{eqnarray}$$

これらが計算できればE-stepでの目的である$p(z_i|x_i)$ として

$$\begin{eqnarray} p(z_i=A|x_i) &=& p(z_i=A,x_i)/p(x_i) \\ p(z_i=B|x_i) &=& p(z_i=B,x_i)/p(x_i) \end{eqnarray}$$

右辺は(2)と(1)から計算できる。

実際にはこれを計算する:

$$\begin{eqnarray} p(z_i=A,x_i) &=& p(x_i|z_i=A)p(z_i=A)\\ p(z_1=A,x_1) &=& p(x_1|z_1=A)p(z_1=A) \end{eqnarray}$$

$p(x_1|z_1=A)$ は $z_1=A$ のときの $p(x_1)=[5\ 5]$ が起こる尤度だから

px1z1A = binopdf(X(1,1), X(1,1)+X(1,2), theta_Aest)
px1z1A =       0.20066

初期設定より、$p(z_1=A) = 0.5$

よって $p(z_1=A,x_1) = p(x_1|z_1=A)p(z_1=A)$ は以下のように計算できる。

pz1Ax1 = px1z1A .* pziAest(1)
pz1Ax1 =       0.10033

同様にして、$p(x_1|z_1=B)$ は $z_1=B$ のときの $p(x_1)=[5\ 5]$ が起こる尤度だから

px1z1B = binopdf(X(1,2), X(1,1)+X(1,2), theta_Best)
px1z1B =       0.24609

初期設定より、$p(z_1=B) = 0.5$

よって $p(z_1=B,x_1) = p(x_1|z_1=B)p(z_1=B)$ は以下のように計算できる。

pz1Bx1 = px1z1B .* pziBest(1)
pz1Bx1 =       0.12305

まとめると $p(z_1,x_1) = [0.10033\ 0.12305]$ となる。

同様にして $x_1,x_2,... x_5$ について計算すると

$$p(z_i=A,x_i) = p(x_i|z_i=A)p(z_i=A)\\ p(z_i=B,x_i) = p(x_i|z_i=B)p(z_i=B)$$

pziAxi = ([binopdf(X(:,1), X(:,1)+X(:,2), theta_Aest)]) .* pziAest;
pziBxi = ([binopdf(X(:,2), X(:,1)+X(:,2), theta_Best)]) .* pziBest;
[pziAxi pziBxi]
ans =
      0.10033      0.12305
     0.020155    0.0048828
     0.060466     0.021973
     0.055738      0.10254
       0.1075     0.058594

$p(x_i)$ を計算するには上で計算した $p(z_i=A,x_i)$ と $p(z_i=B,x_i)$ があればよい。

$$p(xi) = p(z_i=A,x_i) +p(z_i=B,x_i) \tag{1}$$

より

pxi = pziAxi + pziBxi
pxi =
      0.22338
     0.025038
     0.082439
      0.15828
      0.16609

これが$z_i$ の推定値

$$p(zi=A|xi) = p(zi=A,xi) / p(xi)\\ p(zi=B|xi) = p(zi=B,xi) / p(xi)$$

pziAest = pziAxi ./ pxi;
pziBest = pziBxi ./ pxi;
[pziAest pziBest]
ans =
      0.44915      0.55085
      0.80499      0.19501
      0.73347      0.26653
      0.35216      0.64784
      0.64722      0.35278

(3) M-step: $\theta_{Aest}, \theta_{Best}$ の推定

coinAが出る期待値は

$$E(z_i=A|x_i) = p(z_i=A|x_i) p(z_i=A) = p(z_i=A|x_i)p(x_i)$$

X_Aest = X .* pziAest;
X_Best = X .* pziBest;
sumX_Aest = sum(X_Aest,1);
sumX_Best = sum(X_Best,1);
theta_Aest = sumX_Aest(1) / (sumX_Aest(1) + sumX_Aest(2));
theta_Best = sumX_Best(1) / (sumX_Best(1) + sumX_Best(2));
[theta_Aest theta_Best]
ans =       0.71301      0.58134

(4) 収束条件を見てストップするか判断

thetaが収束したかどうかを確認して、E-Mの繰り返しを止めるか決める。

ではE-Mを20回繰り返して、thetaが収束するかどうか見てみることにしよう。

clear all
X = [...
    5 5;
    9 1;
    8 2;
    4 6;
    7 3];
theta_Aest = 0.60;
theta_Best = 0.50;
pziAest = [...
    0.5;
    0.5;
    0.5;
    0.5;
    0.5];
pziBest = 1 - pziAest;
theta_est =[theta_Aest theta_Best];
for ii=1:20
    pziAxi = ([binopdf(X(:,1), X(:,1)+X(:,2), theta_Aest)]) .* pziAest;
    pziBxi = ([binopdf(X(:,2), X(:,1)+X(:,2), theta_Best)]) .* pziBest;
    pxi = pziAxi + pziBxi;
    pziAest = pziAxi ./ pxi;
    pziBest = pziBxi ./ pxi;
    X_Aest = X .* pziAest;
    X_Best = X .* pziBest;
    sumX_Aest = sum(X_Aest,1);
    sumX_Best = sum(X_Best,1);
    theta_Aest = sumX_Aest(1) / (sumX_Aest(1) + sumX_Aest(2));
    theta_Best = sumX_Best(1) / (sumX_Best(1) + sumX_Best(2));
    theta_est = [theta_est; theta_Aest theta_Best];
end
figure;
hold on
plot(0:size(theta_est)-1, theta_est(:,1), 'ro-')
plot(0:size(theta_est)-1, theta_est(:,2), 'bo-')

このように、Zが既知での最尤法のときに計算した値(thetaA = 0.80, thetaB = 0.45)と同じ値に収束していることが分かる。

ところでNature BiotechnologyのPrimerの記事では10回繰り返したところでthetaの推定値が thetaA = 0.80, thetaB = 0.52 となってる。たぶんPrimerの記事のほうがなんか間違えてると思う。

お勧めエントリ

• 細胞外電極はなにを見ているか(1) 20080727 (2) リニューアル版 20081107
• 総説 長期記憶の脳内メカニズム 20100909
• 駒場講義2013 「意識の科学的研究 - 盲視を起点に」20130626
• 駒場講義2012レジメ 意識と注意の脳内メカニズム(1) 注意 20121010 (2) 意識 20121011
• 視覚、注意、言語で3*2の背側、腹側経路説 20140119
• 脳科学辞典の項目書いた 「盲視」 20130407
• 脳科学辞典の項目書いた 「気づき」 20130228
• 脳科学辞典の項目書いた 「サリエンシー」 20121224
• 脳科学辞典の項目書いた 「マイクロサッケード」 20121227
• 盲視でおこる「なにかあるかんじ」 20110126
• DKL色空間についてまとめ 20090113
• 科学基礎論学会 秋の研究例会 ワークショップ「意識の神経科学と神経現象学」レジメ 20131102
• ギャラガー＆ザハヴィ『現象学的な心』合評会レジメ 20130628
• Marrのrepresentationとprocessをベイトソン流に解釈する (1) 20100317 (2) 20100317
• 半側空間無視と同名半盲とは区別できるか？(1) 20080220 (2) 半側空間無視の原因部位は？ 20080221
• MarrのVisionの最初と最後だけを読む 20071213

月別過去ログ

[2024]

5|

4|

3|

2|

1|

[2023]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2022]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2021]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2020]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2019]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2018]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2017]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2016]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2015]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2014]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2013]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2012]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2011]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2010]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2009]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2008]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2007]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2006]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2005]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2004]

12|

11|

10|

9|

8|

7|

6|

5|

4|

3|

2|

1|

[2003]

12|

11|

8|

7|

[2001]

10|

8|

6|

5|

4|

2|

[2000]

12|

9|

8|

7|

6|

3|

2|

1|

[1999]

9|

8|

7|