一、機器學習中的引數估計問題

    在前面的博文中，如“簡單易學的機器學習演算法——Logistic迴歸”中，採用了極大似然函式對其模型中的引數進行估計，簡單來講即對於一系列樣本，Logistic迴歸問題屬於監督型學習問題，樣本中含有訓練的特徵以及標籤，在Logistic迴歸的引數求解中，通過構造樣本屬於類別和類別的概率：

這樣便能得到Logistic迴歸的屬於不同類別的概率函式：

此時，使用極大似然估計便能夠估計出模型中的引數。但是，如果此時的標籤是未知的，稱為隱變數，如無監督的學習問題，典型的如K-Means聚類演算法，此時不能直接通過極大似然估計估計出模型中的引數。

二、EM演算法簡介

    在上述存在隱變數的問題中，不能直接通過極大似然估計求出模型中的引數，EM演算法是一種解決存在隱含變數優化問題的有效方法。EM演算法是期望極大(Expectation Maximization)演算法的簡稱，EM演算法是一種迭代型的演算法，在每一次的迭代過程中，主要分為兩步：即求期望(Expectation)步驟和最大化(Maximization)步驟。

三、EM演算法推導的準備

1、凸函式

    設是定義在實數域上的函式，如果對於任意的實數，都有

那麼是凸函式。若不是單個實數，而是由實陣列成的向量，此時，如果函式的Hesse矩陣是半正定的，即

那麼是凸函式。特別地，如果或者，那麼稱

2、Jensen不等式

    如果函式是凸函式，是隨機變數，那麼

特別地，如果函式是嚴格凸函式，那麼當且僅當

即隨機變數是常量。

(圖片來自參考文章1)

注：若函式是凹函式，上述的符號相反。

3、數學期望

3.1隨機變數的期望

   設離散型隨機變數的概率分佈為：

其中，，如果絕對收斂，則稱為的數學期望，記為，即：

   若連續型隨機變數的概率密度函式為，則數學期望為：

3.2隨機變數函式的數學期望

   設是隨機變數的函式，即，若是離散型隨機變數，概率分佈為：

則：

   若是連續型隨機變數，概率密度函式為，則

四、EM演算法的求解過程

    假設表示觀測變數，

表示潛變數，則此時即為完全資料，的似然函式為，其中，為需要估計的引數，那麼對於完全資料，的似然函式為。

    構建好似然函式，對於給定的觀測資料，為了估計引數，我們可以使用極大似然估計的方法對其進行估計。因為變數是未知的，我們只能對的似然函式為進行極大似然估計，即需要極大化：

上述式子中無法直接對求極大值，因為在函式中存在隱變數，即未知變數。若此時，我們能夠確定隱變數的值，便能夠求出的極大值，可以用過不斷的修改隱變數的值，得到新的的極大值。這便是EM演算法的思路。通過迭代的方式求出引數。

    首先我們需要對引數賦初值，進行迭代運算，假設第次迭代後引數的值為，此時的log似然函式為，即：

表示的是的期望，其中，表示的是隱變數滿足的某種分佈。這樣，上式的值取決於和的概率。在迭代的過程中，調整這兩個概率，使得下界不斷的上升，這樣就能求得的極大值。注意，當等式成立時，說明此時已經等價於。由Jensen不等式可知，等式成立的條件是隨機變數是常數，即：

至此，我們得出了隱變數滿足的分佈的形式。這就是

步。在確定了後，調整引數使得取得極大，這便是

1. 初始化引數，開始迭代；
2. E步：假設為第次迭代引數的估計值，則在第次迭代中，計算：
3. M步：求使極大化的，確定第次的引數的估計值：

五、EM演算法的收斂性保證

迭代的過程能否保證最後找到的就是最大的似然函式值呢？即需要證明在整個迭代的過程中，極大似然估計是單調增加的。假定和是EM演算法的第次和第次迭代後的結果，選定，進行迭代：

1. E步：
2. M步：

固定，將看成變數：

六、利用EM演算法引數求解例項

    假設有有一批資料分別是由兩個正態分佈：

產生，其中，和未知，。但是不知道具體的是第產生，即可以使用和表示。這是一個典型的涉及到隱藏變數的例子，隱藏變數為和。可以使用EM演算法對引數進行估計。

1. 首先是初始化和；
2. E步：，即求資料是由第個分佈產生的概率：
3. M步：，即計算最大的期望值。然而我們要求的引數是均值，可以通過如下的方式估計：

Python程式碼

#coding:UTF-8
'''
Created on 2015年6月7日

@author: zhaozhiyong
'''
from __future__ import division
from numpy import *
import math as mt
#首先生成一些用於測試的樣本
#指定兩個高斯分佈的引數，這兩個高斯分佈的方差相同
sigma = 6
miu_1 = 40
miu_2 = 20

#隨機均勻選擇兩個高斯分佈，用於生成樣本值
N = 1000
X = zeros((1, N))
for i in xrange(N):
    if random.random() > 0.5:#使用的是numpy模組中的random
        X[0, i] = random.randn() * sigma + miu_1
    else:
        X[0, i] = random.randn() * sigma + miu_2

#上述步驟已經生成樣本
#對生成的樣本，使用EM演算法計算其均值miu

#取miu的初始值
k = 2
miu = random.random((1, k))
#miu = mat([40.0, 20.0])
Expectations = zeros((N, k))

for step in xrange(1000):#設定迭代次數
    #步驟1，計算期望
    for i in xrange(N):
        #計算分母
        denominator = 0
        for j in xrange(k):
            denominator = denominator + mt.exp(-1 / (2 * sigma ** 2) * (X[0, i] - miu[0, j]) ** 2)
        
        #計算分子
        for j in xrange(k):
            numerator = mt.exp(-1 / (2 * sigma ** 2) * (X[0, i] - miu[0, j]) ** 2)
            Expectations[i, j] = numerator / denominator
    
    #步驟2，求期望的最大
    #oldMiu = miu
    oldMiu = zeros((1, k))
    for j in xrange(k):
        oldMiu[0, j] = miu[0, j]
        numerator = 0
        denominator = 0
        for i in xrange(N):
            numerator = numerator + Expectations[i, j] * X[0, i]
            denominator = denominator + Expectations[i, j]
        miu[0, j] = numerator / denominator
        
    
    #判斷是否滿足要求
    epsilon = 0.0001
    if sum(abs(miu - oldMiu)) < epsilon:
        break
    
    print step
    print miu
    
print miu

[[ 40.49487592  19.96497512]]

參考文章：

1、(EM演算法)The EM Algorithm (http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html)

2、數學期望(http://wenku.baidu.com/view/915a9c1ec5da50e2524d7f08.html?re=view)