EM演算法的個人理解和實驗

引言

很早以前就接觸了EM演算法，但是發覺自己好像有理解不正確的地方，於是自己推導了一遍。才疏學淺，有錯誤之處還請各位大佬指正。

EM演算法的通俗解釋

Expectation Maximization 的演算法的本質，就是對一組數據 ${x_1, x_2, x_3, ...}$ ，假設他們滿足某種分佈 , 求出這個最優的 , 使得的值最大。在不存在隱變數的情況下，也就是說，，則這個問題就是個簡單的求函數極值的問題，也就是最普通的參數估計問題，相信大家概率論的課上都學過了。

但是在實際情況下，除去分佈的參數，還有一些隱變數，也決定了的分佈。一個簡單的例子就是對於很多數據點 , 其中的任一個點都是兩個正態分佈中的其中一個。那麼每一個屬於哪一個正態分佈，產生的標籤 $c_i in {0, 1}$ 就是隱變數。

但是我們是不是也可以把隱變數當做分佈的參數呢？也就是說從兩個正態分佈裡面取樣產生的新的分佈的參數包含 ${ mu_1, sigma_1, mu_2, sigma_2, c_1, c_2, c_3, ....}?$ ，搜索整個參數空間，理論上也可以找到一組最佳的參數，使得最大。

但是這樣就和我們的目標有所違背：我們的目標是找到的最大值，並不是的最大值。而且一般的函數也十分複雜，很可能無法求出極值。

通過貝葉斯法則，我們有： $p(x| heta) = dfrac{p(x, h| heta)}{p(h|x, heta)}?$ ，但是這個值和有關，所以很難直接對這個公式下手。EM演算法提供了一種巧妙的方法：

對上面的等式求對數，得到: , 同時又有 , 因此把這個對 $h?$ 加權求和的算符作用在左右兩邊，很顯然有 , 出現了很熟悉的的形式！貌似可以利用資訊理論裡面的那些不等式求極值了！

記住左邊求和項，對於任何都成立。因為等式左邊不含 , 於是我們可以把上式改成。觀察等號左邊的第二項，就是在條件下的信息熵。資訊理論告訴我們信息熵是最短編碼長度。那麼如果變數的分佈發生改變，但編碼方式不變（也就是說左邊的變化了），必然會增加平均碼長。也就是說。當然其實這本質就是「相對熵必然為正」的 Gibbs不等式。

這告訴我們，把右邊的換成，有，等號右邊的第二項是一定增加的，所以我們只要保證增加就行。

這就是EM演算法的精髓，通過對於當前的，找到一個，使得取到最大值，不斷迭代，最終找到最優的。

一個簡單的例子

就用上文的兩個正態分佈的例子，對於數據點 ${x_1, x_2, ... ,x_n}?$ ，他們屬於兩個正態分佈之一，參數分別為 ${mu_1, sigma_1, mu_2, sigma_2}?$ ，數據的標籤為 ${c_1,c_2,...,c_n}?$ 。表示表示。這時候我們把數據的標籤作為隱變數。

假設目前已經進行到 t 步，即當前的參數分別為 ${mu_{1,t}, sigma_{1,t}, mu_{2,t}, sigma_{2,t}}?$ , 此時有我們要求出優化目標求和號裡面的左邊部分 $p(h|x, heta_t) = dfrac{p(h,x, heta_t)}{p(x, heta_t)} = dfrac{p(h, x| heta_t)}{p(x| heta_t)}?$ 。

這個式子依然難以計算，但是注意到

，左邊的一般都是有函數表達式的，右邊的，在這裡我們認為隱變數 $h$ 和參數是無關的，並且認為的分佈由另一組參數決定，即：。在這裡，我們很容易認為單個數據的標籤滿足伯努利分佈，即：有的概率標籤為 1，表示該數據點屬於第一個正態分佈；反之則屬於第二個。所以也要加入參數中，參數集合變為： ${mu_1, sigma_1, mu_2, sigma_2, au}$

於是我們有：，

$sumlimits_h p(h|x, heta_t) log p(x, h| heta) = sumlimits_h dfrac{p(h)cdot p(x|h, heta_t)}{sumlimits_h p(h) cdot p(x|h, heta_t)} log p(x|h, heta)cdot p(h)?$

對於任何一個數據點，對於隱變數 , 只需要考慮的值，因為別的數據點處於哪一個分佈對沒有影響。因此可以把隱變數的集合拆成 $H_1 = {h|h_i = 1}, H_2 = {h|h_i = 2}$ 兩個集合。這樣我們就可以把公式拆成：

$dfrac{1}{sumlimits_{hin H_1} p(h) p(x_i|h_i = 1, heta_t) +sumlimits_{hin H_2}{p(h)p(x_i|h_i=2, heta_t)}} cdot (sumlimits_{h in H_1} p(h) p(x_i|h_i = 1, heta_t) log p(h) p(x_i|h_i = 1, heta) + sumlimits_{hin H_2} p(h) p(x_i|h_i = 2, heta_t) log p(h) p(x_i|h_i = 2, heta))$

顯然有： $sumlimits_{hin H_1} p(h) = au, sumlimits_{hin H_2} p(h) = 1- au?$

最後我們需要考慮的是，上面的討論僅僅是最大化某一個樣本的 $log p(x| heta)?$ 的值的，考慮到我們有大量的樣本，實際的優化目標 $p(X| heta) = prodlimits_{x in X} p(x| heta) ?$ ，取對數後相加即可。因此在上式的左側還要加上 $sumlimits_{x_i in X}?$

於是最終的式子為：

$sumlimits_{x_i in X} sumlimits_{h_i in {1,2} } p(h_i|x_i, heta_t) log au_{h_i} p(x_i|h_i, heta)?$

其中， $au_1 = au, au_2 = 1 - au, p(h_i|x_i, heta_t) = dfrac{ au_{h_i}cdot p(x_i|h_i, heta_t)}{sumlimits_{h_j in {1, 2}} au_{h_j}cdot p(x_i|h_j, heta_t)}$

考慮到上式的右端 $log au_{h_i}p(x_i|h_i, heta)?$ 中的 $au_{h_i}?$ 與無關，因此可以忽略之。再看這個式子，是不是很像

$sumlimits_{x in X} w_x log p(x| heta)?$ 的形式？因為這是對數形式，我們還原出其原始式子： $prodlimits_{x in X} p(x| heta)^{w_x}?$ 很顯然這就是一個概率值，可以理解為樣本出現了次數的概率。注意到不同的對應的屬於不同的正態分佈，所以我們要把關於的那個求和的兩項分別拿出來求極值。

拿出這一項， $p(x_i|h_i, heta_t) = dfrac{1}{sqrt{2pi}sigma_1} e^{-dfrac{(x_i - mu_1)}{2sigma_1^2}} = f(x_i;mu_1,sigma_1)?$ 。因此這就是個正態分佈的加權最大似然估計（Weighted MLE for Normal Distribution），很容易可以得到：

$mu_{1,t} = dfrac{sumlimits_idfrac{ au_1 cdot f(x_i;mu_1, sigma_1)}{ au_1 cdot f(x_i;mu_1, sigma_1) + au_2 cdot f(x_i;mu_2, sigma_2)} cdot x_i}{sumlimits_idfrac{ au_1 cdot f(x_i;mu_1, sigma_1)}{ au_1 cdot f(x_i;mu_1, sigma_1) + au_2}} = dfrac{sumlimits_i T_{i,1} x_i}{sumlimits_i T_{i, 1}}?$

其中， $T_{i,j} = dfrac{ au_j cdot f(x_i;mu_j, sigma_j)}{ au_1 cdot f(x_i;mu_1, sigma_1) + au_2 cdot f(x_i;mu_2, sigma_2)}quad j in {1,2}?$

$sigma_{1,t}^2 = dfrac{sumlimits_i T_{i,1}(x_i - mu_{1,t})^2}{sumlimits_i T_{i, 1}}$

當然同時，我們還要注意到也是一個參數，因此也需要對此進行更新。此時可以把 $p(h_i|x_i, heta_t) log au_{h_i} p(x_i|h_i, heta)$ 中的裏的看做常數。此時形式變成了

或者是 $au_1^{w_1} (1- au_1)^{w_2}$ 很顯然這就是一個二項分佈的概率表達式，相當於從一堆紅球和黑球中選了若干次，總共有次紅球，次黑球，因此其最大似然估計是 $au_1 = dfrac{w_1}{w_1 + w_2},quad au_2 = dfrac{w_2}{w_1 + w_2}$

這裡 $w_1 = sumlimits_i T_{i,1}$ ，而且可以注意到 $sumlimits_i T_{i,1} + sumlimits_i T_{i,2} = 1$

所以有 $au_1 = sumlimits_i T_{i, 1}$

一些理解

從上面的例子中，我們可以看到某個數據點屬於哪一個分佈，這是一組隱變數，且這組隱變數中的每一個值滿足伯努利分佈。隱變數自身也包含了參數。因此我們可以畫出如下的關係圖：

，，而且是我們能觀測到的數據。EM演算法做的事情就是在把表示成關於的函數時，找出其最大值。

一個簡單的實驗

import numpy from matplotlib import pyplot def em_two_random_normal(mean_1, var_1, mean_2, var_2, size, n_steps = 10): x1s = numpy.random.normal(mean_1, var_1, [size, 2]) x2s = numpy.random.normal(mean_2, var_2, [size, 2]) xs = numpy.vstack([x1s, x2s]) mu_1_s = mean_1 + numpy.random.uniform(0, 8, [1, 2]) mu_2_s = mean_2 + numpy.random.uniform(-8, 0, [1, 2]) var_1_s = 1 var_2_s = 1 p_z1 = 0.5 p_z2 = 0.5 mus = [[mu_1_s], [mu_2_s]] vars = [[var_1_s], [var_2_s]] for _ in range(n_steps): p_x_at_para_z1 = numpy.exp(-numpy.sum(numpy.square(xs - mu_1_s), axis=1, keepdims=True) / (2 * var_1_s)) / numpy.sqrt(2 * numpy.pi * var_1_s) p_x_at_para_z2 = numpy.exp( - numpy.sum(numpy.square(xs - mu_2_s), axis=1, keepdims=True) / (2 * var_2_s)) / numpy.sqrt(2 * numpy.pi * var_2_s) p_z1_under_para_x = p_z1 * p_x_at_para_z1 / (p_z1 * p_x_at_para_z1 + p_z2 * p_x_at_para_z2) p_z2_under_para_x = p_z2 * p_x_at_para_z2 / (p_z1 * p_x_at_para_z1 + p_z2 * p_x_at_para_z2) p_z1 = numpy.mean(p_z1_under_para_x) p_z2 = numpy.mean(p_z2_under_para_x) mu_1_s = numpy.sum(p_z1_under_para_x * xs, axis=0, keepdims=True) / numpy.sum(p_z1_under_para_x) mu_2_s = numpy.sum(p_z2_under_para_x * xs, axis=0, keepdims=True) / numpy.sum(p_z2_under_para_x) var_1_s = numpy.sum(p_z1_under_para_x * numpy.sum(numpy.square((xs - mu_1_s)), axis=1, keepdims=True)) / numpy.sum(p_z1_under_para_x) var_2_s = numpy.sum(p_z2_under_para_x * numpy.sum(numpy.square((xs - mu_2_s)), axis=1, keepdims=True)) / numpy.sum(p_z2_under_para_x) mus[0].append(mu_1_s) mus[1].append(mu_2_s) vars[0].append(var_1_s) vars[1].append(var_2_s)

pyplot.scatter(xs[:, 0], xs[:, 1])
for i in range(n_steps):
pyplot.gca().add_artist(pyplot.Circle([mus[0][i][0, 0], mus[0][i][0, 1]], numpy.sqrt(vars[0][i]), color=red, fill=False))
pyplot.gca().add_artist(pyplot.Circle([mus[1][i][0, 0], mus[1][i][0, 1]], numpy.sqrt(vars[1][i]), color=blue, fill=False))
pyplot.xlim(0, 20)
pyplot.ylim(0, 20)
pyplot.show()