EM Algorithm

時間 2019-12-07

標籤 algorithm 欄目 CSS 简体版

原文原文鏈接

#Expectation Maximization Algorithm EM算法和以前學的都不太同樣，EM算法更多的是一種思想，因此後面用幾個例子講解，同時也會重點講解GMM高斯混合模型。 ###①極大似然估計極大似然估計這裏面用的比較多。假設咱們想要知道咱們學生身高的分佈，首先先假設這些學生都是符合高斯分佈$$N(μ, σ^2)$$咱們要作的就是要估計這兩個參數究竟是多少。學生這麼多，挨個挨個來確定是不切實際的，因此天然就是抽樣了。爲了統計學生身高，咱們抽樣200我的組成樣本git

咱們須要估計的參數$$θ = [μ, σ]^T$$ 首先估計一下抽到這兩百人的機率一共是多少，抽到男生A的機率github

P(x_A|θ)$$抽到學生B的機率$$P(x_B|θ)$$因此同時抽到這兩個學生的機率就是$$P(x_A|θ) * P(x_B|θ)

那麼同時抽到這200個學生的G機率$$L(θ) = L(x_1,x_2,x_3...,x_200;θ) = \prod_{i=1}^{200}P(x_i|θ)$$ 最後再取一個對數就行了：算法

####notation和log 上面有一條公式裏面是同時存在了;和|，這兩個符號差異其實有點大的。**|通常咱們是用來表示條件機率，好比$$P(x|θ)$$就是表示x在θ的條件下發生的機率。框架

也是一個意思。分號;表示的就是表示後面的是待估計的參數，也就是說P(x;θ)意思就是後面的θ是須要估計的參數而不是條件，因此|也有另外一層意思，若是不是表示條件機率，那麼就是表示後面有待估計參數。固然是在|不表示條件機率的狀況下。** 這兩種表示法是源於兩種學派的不一樣理解：頻率派認爲參數爲固定的值，是指真實世界中，參數值就是某個定值。貝葉斯派認爲參數是隨機變量，是指取這個值是有必定機率的。固然，不管是;仍是|，他們都是表示條件機率的意思，只不過這兩個學派理解不同而已。 notation的問題解決了以後就是log的問題了，爲何須要log化，講道理，是不須要的。可是求log有這麼幾個好處： 1.計算簡單，累乘是很難計算的，log以後能夠變換成累加。 2.機率累乘是會出現數字很是小的狀況，log以後就能夠避免這種狀況。 3.log以後函數的梯度方向是沒有變化的，對於函數優化的方向影響很小。dom

似然函數的執行步驟： 1.獲得似然函數 2.取對數整理 3.求導數，另導數爲零 4.解方程獲得解函數

###②Jensen不等式首先引出凸函數的概念$$f(x)^{''} > 0$$那麼就是凸函數，因此它的圖像就是一個勾形的，看起來是一個凹函數，其實是凸函數。學習

凸函數的性質很顯而易見了$$E[f(x)] >= f(E(x))$$

其實很明顯的，看圖就知道了，E(x)其實就是a到b的平均值，上面的結論很容易證明。那麼若是想要取到等號，須要改變上面？取等號的意思就是相切，相切的意思就是a = b，既然a = b了，那麼x天然就是常數了，因此當且僅當$$P(x = E(x)) == 1$$若是是凹函數，那麼就是方向相反了。 ###③EM算法的推導正常來看先是要引入一個最大似然函數：$$L(θ) = \sum_{i=1}^{n}logP(x;θ)$$但這樣實際上是和難求的，P(x|θ)徹底混在了一塊兒，根本求不出來，因此咱們要引入一個輔助變量z。

####隱變量Z 隱變量是觀測不到的，好比作一個抽樣，從3個袋子裏面抽取小球。而抽取這些小球的過程是不可見的，抽取的過程其實就是隱變量，而這些隱變量，也就是過程能夠告訴咱們這個x是從哪一個袋子來的，由此來區分。這個隱變量和HMM裏面的隱含序列不是一個東西，隱含序列是當前x的狀態，而不是抽取x的過程。因此在這裏，通俗點講，這個z就是用來找到x的組類的，也就是說z來告訴這個x你是屬於哪一組的。另外須要注意的是，隱變量是不能夠隨便添加的，添加隱變量以後不能影響邊緣機率。也就是說，原來是P(x)，添加以後就是P(x,z)，那麼必須有： $P(x) = \int_zP(x, z){\rm dz}$ 優化

因此咱們引入隱變量的緣由是爲了轉化成和這幾個高斯模型相關的式子，不然無從下手。化簡一下上式子：$$L(θ) = \sum_{i=1}^{n}logP(x;θ)=\sum_{i=1}^{n}log\sum_zP(x, z;θ)$$既然z能夠指定x，那麼咱們只須要求解出z就行了。注意上面凸函數所提到的一個指望性質，這裏就可使用了。由於雖然優化了上面的式子，仍是不能求出來，由於z變量實在是太抽象了，找不到一個合適的公式來表示它。EM的一個方法就是用優化下界函數的方法來達到優化目標函數的目的。既然z很抽象，那麼咱們就須要一個轉變一下。對於每個樣例x都會對應一個z，那麼假設一個分佈Q(z)是知足了z的分佈的，而Q(z)知足的條件是 $\sum_zQ_i(z) == 1,Q_i(z) > 0$ Qi意味着每個x對應的z都會對應着一個Q了，這裏有點複雜，再詳細解釋一下。一個x對應一組z，z是一個向量，可是每個z又會分別對應一個一個分佈Q。覺得最後獲得的z不會是一個數字，而是一個機率，也就是說Q(z)獲得的是這個x樣例屬於這個類別的機率是多少。而z的數量，一個是當前有多少個分佈混合在一塊兒的數量。再梳理一下：如今的樣本是xi，那麼每個xi將會對應着一組的z，每個xi同時也會對應着一個分佈Qi，z其實就是反應了這個樣本是來自於哪一個分佈的。好比這個x是A1分佈作了3，A2分佈作了5，那麼z可能就是={3,5}。因此Qi(z)獲得的是這個x屬於這些個分佈的機率，也就是說這些分佈對x作了多少百分比的功，天然就是要等於1了。還要注意的是，上面的 $\sum_zQ_i(z) = 1$ 這個並不能獲得Qi(z)就是分佈對x作了多少功的結論，獲得這個結論是後面下界函數與目標函數相等獲得的。這裏只是知道了總和等於1，由於是分佈的總和嘛。如今就到了公式的化簡：$$\sum_ilogP(x^{(i)};θ) = \sum_ilog\sum_{z^{(i)}}P(x^{(i)},z^{i};θ) = \sum_ilog\sum_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})}$$ 仔細看一下這個式子$$\sum_zQ_i(z^{(i)})\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})}$$這個式子其實就是求 $\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})}$ 的指望，假設 $\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})} = x$ ，那麼能夠利用上面。因而化簡：.net

\sum_ilog\sum_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})} >= \sum_i\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})}$$這個時候就獲得了下界函數，上面也講過了，想要相等，天然就是x要是常數，因此$$\frac{P(x^{i},z^{i};θ)}{Q_i(z^{(i)})} = c$$既然$\sum_zQ_i(z) = 1$，並且z也是同樣的，由於一個樣本嘛。因此上下加和（若是是離散的，那就sum一下，連續的那就積分，這裏是離散的，因此就是sum一下）。因而有$$\frac{\sum_zP(x^i,z^i;θ)}{\sum_zQ_i(z^i)} = c

因而有： 3d

因此Q(z)計算的是後驗機率。計算在當前參數θ和已經抽樣到x的條件下，這個x是從z來的機率。其實就是z對x作了多少貢獻。因此整個EM算法步驟就很清晰了：

####EM算法計算步驟： E-step：對於每個，求 $Q_i^{(t)}(z) = P(z^i|x^i;θ^{(t)})$ M-step： $θ = argmax_θ\sum_i\sum_{z^{(i)}}Q_i(z^i)log\frac{P(x^i,z^i;θ)}{Q_i(z^i)}$ 這時候就可使用求導迭代的方法求解了。

這就是整一個EM算法的框架了，能夠看到其實沒有比較具體的算法，大體上就是一個框架。那麼問題來了，怎麼樣證實這東西是一個收斂的？？

####證實EM算法的收斂性既然是極大似然估計，那麼須要證實的天然就是 $L(θ^{t}) <= L(θ^{(t+1)})$ ，也就是說極大似然估計要單調遞增。由於每一次都是極大，那麼隨着時間增長，天然就是要增大的了。當給定一個 $θ^{(t)}$ 的時候，至關於t時間的EM算法進行了E步了，因此有：

L(θ^{t}) = \sum_i\sum_{z^{(i)}}Q_i^{t}(z^i)log\frac{P(x^i,z^i;θ^{t})}{Q_i(z^i)}$$而後就走M步，M步以前尚未求最大，因此:$$L(θ^{(t+1)}) >= \sum_i\sum_{z^i}Q_i^{(t)}(z^i)log\frac{P(x^i,z^i;θ^{(t+1)})}{Q_i(z^i)}$$**這一步有點複雜，這裏是M-step了，$θ^{(t+1)}$是L(θ)求極大值獲得的，而最重要的是L和後面下界函數的關係。這兩個東西是不相等的，下界函數等於目標函數的條件是x爲constant，x == constant意味着Q得求出來，這裏看這個Q貌似是出來了，可是這個Q是$Q^t_i$，是上一個時刻的，而不是這個時刻的，因此t+1時刻Q尚未被固定住，也就是沒有知足x == constant的條件，因此下界函數小於等於目標函數。接下來就簡單了，直接把$θ^{(t+1)}$換成$θ^{(t)}$就行了。**![](https://upload-images.jianshu.io/upload_images/10624272-e4c4d749e2aa0ded.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
這樣就證實了收斂性。

>>####EM algorithm的優化方法：
>>以前討論過，這種方法是迭代，使用極大似然估計和迭代的方法來進行優化，但實際上這更像是一種座標上升的方法：
![](https://upload-images.jianshu.io/upload_images/10624272-544bd237ec99baf9.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
一次固定一個變量，對另外的求極值，最後逐步逼近極值。對應到EM上，E步：固定 θ，優化Q；M步：固定 Q，優化 θ；交替將極值推向極大。Kmeans也是這樣更新的，固定中心點，優化Ein；優化完成，更新中心點。SMO也是，固定$α_{3-n}$更新α1和α2。

#④GMM的推導
能夠直接把高斯混合模型代入EM框架裏面。
**存在多個高斯分佈混合生成了一堆數據X，取各個高斯分佈的機率是$φ_1,φ_2,φ_3,...$**，第i個高斯分佈的均值是$μ_i$，方差是$σ$，求法φ，μ，σ。
按照套路，第一個E-step求出Q，因而有：$$w_j^{(i)} = Q_i(z^{(i)} = j) = P(z^{i} = j |x^i;φ,μ,Σ)

意思就是求出第i個樣本屬於第j個分佈的機率是多少。以後就是M-step了，就是化簡了：

\sum_{i=1}^m\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^i,z^i;φ,μ,Σ)}{Q_i(z^{(i)})}=\sum_{i=1}^m\sum_{j=1}^kQ_i(z^i = j)log\frac{P(x^i|z^i = j;μ,Σ)P(z^i = j;φ)}{Q_i(z^i = j)}

這裏可能須要解釋一下，根據至於條件，由於很明顯，z是隱變量，只是指明瞭x是屬於哪一個類別，和μ，Σ沒有什麼關係，因此直接忽略那兩個參數了，因此P(z)是沒有那兩個參數的，z是表明了分佈，因此每個分佈的機率確定是包括了，因此就只有一個機率的參數。P(x|z)是自己的機率，就是已經知道分佈是那個了，求屬於這個分佈的機率是多少，既然已經選定了分佈那麼天然就不須要再看φ了，由於φ是各個分佈的機率。

以後就是求導數了：

導數等於0，因而有：$$μ_j = \frac{\sum_{i=1}^mw_j^{(i)}x^{(i)}}{\sum_{i=1}^{m}w_j^{(i)}}$$μ就出來了。對於φ，那就簡單了，實際上須要優化的式子：

\sum_{i=1}^{m}\sum_{j=1}^{k}w_j^{i}logφ_j$$求導等於0。
要注意的是φ還有一個條件：$\sum_{j=1}^{k}φ_j = 1$。有條件，那麼問題來了，不能直接求導，因此只能用拉格朗日乘子法求解，因而構造式子：$L(φ) = \sum_{i=1}^{m}\sum_{j=1}^{k}w_j^ilogφ_j + β(\sum_{j=1}^{k}φ_j - 1)$
接着求導$\sum_{i=1}^{m}\frac{w_j^{(i)}}{φ_j}+β = 0$。因而獲得$φ_j = \frac{\sum_{i=1}^{m}w_j^{i}}{-β}$，再改變一下，兩邊加和，獲得$\sum_{j=1}^{k}φ_j = \frac{\sum_{i=1}^{m}\sum_{j=1}^{k}w_j^{(i)}}{-β}$，這個就簡單了，$\sum_{j=1}^{k}w_j^i = 1$，因此-β = m。那麼：
![](https://upload-images.jianshu.io/upload_images/10624272-e76886d2bd920041.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
Σ也是同樣：
![](https://upload-images.jianshu.io/upload_images/10624272-735036c5f3e898e1.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
這樣整個推導就算結束了。
#⑤硬幣模型
如今有兩個硬幣AB，進行5次試驗每一次投10次，並不知道是哪一個硬幣投的，求兩種硬幣的正面的機率。
首先E-step：
首先先初始化一下，$θ_A = 0.6,θ_B = 0.5$
第一個試驗選中A的機率：
$$P(z = A|x_1;θ) = \frac{P(z = A,x_1|θ)}{P(z = A,x_1|θ)+P(z = B,x_1|θ)} = \frac{(0.6)^5*(0.4)^5}{(0.6)^5*(0.4)^5 + (0.5)^{10}} = 0.45$$一樣求得$P(z = B|x_1;θ) = 0.55$
計算機出每個試驗的機率而後相加求均值。
以後就是M-step了：
![](https://upload-images.jianshu.io/upload_images/10624272-2d21573def70f1f6.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
$μ_j$表示選擇了A的機率，$1-μ_j$表示選擇了B的機率。以後就是求導更新了。
#⑥代碼實現
方差的求解就不玩了，主要就是迭代求解μ和φ的值了。
首先是生成數據，4個高斯分佈，每個高斯分佈的sigma都是同樣的，不同的只有μ和α，也就是φ，習慣上把前面的一個參數叫作權值，因此用α來表示。
```
def generate_data(sigma, N, mu1, mu2, mu3, mu4, alpha):
global X
X = np.zeros((N, 2))
X = np.matrix(X)
global mu
mu = np.random.random((4, 2))
mu = np.matrix(mu)
global expect
expect = np.zeros((N, 4))
global alphas
alphas = [0.25, 0.25, 0.25, 0.25]
for i in range(N):
if np.random.random(1) < 0.1:
X[i, :] = np.random.multivariate_normal(mu1, sigma, 1)
elif 0.1 <= np.random.random(1) < 0.3:
X[i, :] = np.random.multivariate_normal(mu2, sigma, 1)
elif 0.3 <= np.random.random(1) < 0.6:
X[i, :] = np.random.multivariate_normal(mu3, sigma, 1)
else:
X[i, :] = np.random.multivariate_normal(mu4, sigma, 1)
plt.title('Generator')
plt.scatter(X[:, 0].tolist(), X[:, 1].tolist(), c = 'b')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
```
這四個模型的比例分別是1:2:3:4，使用EM來找到他們屬於的類別。
![](https://upload-images.jianshu.io/upload_images/10624272-ed6bdb6587db87d9.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
其實若是用kmeans聚類的話更加快速，可是這裏仍是用EM。
E-step：
```
def e_step(sigma, k, N):
global X
global mu
global expect
global alphas
for i in range(N):
W = 0
for j in range(k):
W += alphas[j] * math.exp(-(X[i, :] - mu[j, :]) * sigma.I * np.transpose(X[i, :] - mu[j, :])) / np.sqrt(np.linalg.det(sigma))
for j in range(k):
w = math.exp(-(X[i, :] - mu[j, :]) * sigma.I * np.transpose(X[i, :] - mu[j, :])) / np.sqrt(np.linalg.det(sigma))
expect[i, j] = alphas[j]*w/W
pass
```
就是按照公式來求解w便可，求解每個分佈對樣本點作了多少的功，以後求單個樣本點求比例。
M-step：
```
def m_step(k, N):
global expect
global X
global alphas
for j in range(k):
mathor = 0
son = 0
for i in range(N):
son += expect[i, j]*X[i, :]
mathor += expect[i, j]
mu[j, :] = son / mathor
alphas[j] = mathor / N
```
直接按照公式優化便可。
```
if __name__ == '__main__':
iterNum = 1000
N = 500
k = 4
probility = np.zeros(N)
mu1 = [5, 35]
mu2 = [30, 40]
mu3 = [20, 20]
mu4 = [45, 15]
sigma = np.matrix([[30, 0], [0, 30]])
alpha = [0.1, 0.2, 0.3, 0.4]
generate_data(sigma, N, mu1, mu2, mu3, mu4, alpha)
for i in range(iterNum):
print('iterNum : ', i)
err = 0
err_alpha = 0
Old_mu = copy.deepcopy(mu)
Old_alpha = copy.deepcopy(alphas)
e_step(sigma, k, N)
m_step(k, N)
for z in range(k):
err += (abs(Old_mu[z, 0] - mu[z, 0]) + abs(Old_mu[z, 1] - mu[z, 1]))
err_alpha += abs(Old_alpha[z] - alphas[z])
if (err <= 0.001) and (err_alpha < 0.001):
print(err, err_alpha)
break
color = ['blue', 'red', 'yellow', 'green']
markers = ['<', 'x', 'o', '>']
order = np.zeros(N)
for i in range(N):
for j in range(k):
if expect[i, j] == max(expect[i, ]):
order[i] = j
plt.scatter(X[i, 0], X[i, 1], c = color[int(order[i])], alpha=0.5, marker=markers[int(order[i])])
plt.show()
print('standedμ:',mu4, mu3, mu2, mu1)
print('standedα:',alpha)
print('new μ:', mu)
print('new α:',alphas)
```
運行函數。看看結果：
![](https://upload-images.jianshu.io/upload_images/10624272-a85506d016bce0bd.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
![](https://upload-images.jianshu.io/upload_images/10624272-5cf9e38f6ac2a4e6.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
結果其實仍是相差不大。達到預期。

#⑦EM的另外一種理解方法
上面所講的其實只是一種理解方法，在李航老師的統計學習方法裏面是另外一種比較厲害的解法：
**1.E-step：求出Q函數。
2.M-step：利用Q函數求極大值。**
**其實這兩種方法是徹底同樣的，Q函數就是下界函數，![](https://upload-images.jianshu.io/upload_images/10624272-1c18b5cacb8882f9.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)這個公式和上面的：![](https://upload-images.jianshu.io/upload_images/10624272-94b9cd251a48dcf1.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)是同樣的，至於爲何有下面這個$Q_i(z^i)$![](https://upload-images.jianshu.io/upload_images/10624272-3c3ac6b0ff7f9663.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)這是李航老師書上的，既然θ是已經知道了，那麼天然就能夠去掉，由於在log下面作分母就至關因而常數了，拆開能夠變減號嘛，在前面的就不能夠。回來看看咱們剛剛推導的，實際上是同樣的，下面的那個Q確實能夠去掉，由於是事先知道的了。在使用Jensen不等式的時候，須要假設隱變量服從某種形式的機率分佈，才能夠將推導過程的一部分當作是指望的表達形式從而應用Jensen不等式。然而這個分佈不是隨便指定的。咱們令Jensen不等式取等號的時候，能夠計算出這個分佈其實就是：已知觀測數據的隱變量的後驗機率分佈。因爲求Q函數須要先求出隱變量的後驗機率的指望，所以，這就能夠解釋爲何EM算法的「通俗」理解角度的E步驟是求隱變量的指望了。有時候在用EM算法解決某個具體問題的時候，會發現M步驟極大化的竟然是徹底數據的對數似然函數。這是由於，Q函數雖然是徹底數據的對數似然函數的某種指望，可是求這個指望的過程有時其實就是將隱變量的後驗機率的指望代入就能夠了。所以，本質上咱們其實仍是在求Q函數的極大。**
**事實上，李航老師的Estep求出Q函數，其實就是通俗理解裏面的求出Q(z)並代入下界函數的過程。由於求出Q(z)就至關於這個Q(z)被固定了，能夠去掉了，因而就和李航老師的式子同樣了。有機會再補一些李航老師的推導，勃大精深。**
#⑧總結
**EM和Kmeans算法其實很相似，事實上步驟基本能夠用EM框架來替換，可是Kmeans算法是硬分類，說一是一，可是EM算法不太同樣，是軟分類，百分之幾是那個，百分之幾是這個。**

**缺點也仍是有的：初值敏感，局部最優。由於存在了隱變量，因此致使了直接對x作極大似然是不可行的，log已經在sum的外面了。因此EM算法就轉向了下界函數，而這種方法原本就不保證找到局部最優解。**

**若是將樣本看做觀察值，潛在類別看做是隱藏變量，那麼聚類問題也就是參數估計問題。若是一個目標函數存在多個變量，那麼梯度降低牛頓法這些逼近方法就用不了了。但咱們可使用座標上升方法，固定一個變量，對另一個求導數，而後替換最後逐步逼近極值點。對應到EM算法也是同樣，E步求隱含的z變量，Mstep求解其餘參數。**

>>####最後符上GitHub代碼：
>>**https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/ExpectationMaximization**

相關標籤/搜索

algorithm

algorithm&datastructure

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。