EM算法計(jì)算步驟：

E-step：
對于每一個者祖，求
M-step：

這時候就可以使用求導(dǎo)迭代的方法求解了。

這就是整一個EM算法的框架了皮仁，可以看到其實(shí)沒有比較具體的算法谣蠢，大致上就是一個框架册烈。那么問題來了，怎么樣證明這東西是一個收斂的巨坊？究珊？

證明EM算法的收斂性

既然是極大似然估計(jì)瞬浓，那么需要證明的自然就是宪躯，也就是說極大似然估計(jì)要單調(diào)遞增精置。因?yàn)槊恳淮味际菢O大火欧，那么隨著時間增加霎俩，自然就是要增大的了。
當(dāng)給定一個的時候，相當(dāng)于t時間的EM算法進(jìn)行了E步了贷揽，所以有：
然后就走M(jìn)步，M步之前還沒有求最大，所以:這一步有點(diǎn)復(fù)雜回还，這里是M-step了闻葵，是L(θ)求極大值得到的憎夷，而最重要的是L和后面下界函數(shù)的關(guān)系饮焦。這兩個東西是不相等的，下界函數(shù)等于目標(biāo)函數(shù)的條件是x為constant毯炮，x == constant意味著Q得求出來蛋叼，這里看這個Q貌似是出來了，但是這個Q是，是上一個時刻的匣吊，而不是這個時刻的淀歇，所以t+1時刻Q還沒有被固定住掸哑，也就是沒有滿足x == constant的條件琢岩，所以下界函數(shù)小于等于目標(biāo)函數(shù)拼坎。接下來就簡單了赂摆，直接把換成就好了逊朽。

這樣就證明了收斂性渐逃。

EM algorithm的優(yōu)化方法：

之前討論過舶衬，這種方法是迭代朵栖，使用極大似然估計(jì)和迭代的方法來進(jìn)行優(yōu)化坑质，但實(shí)際上這更像是一種坐標(biāo)上升的方法：

一次固定一個變量驻仅，對另外的求極值昧碉，最后逐步逼近極值今妄。對應(yīng)到EM上受神，E步：固定 θ排吴，優(yōu)化Q；M步：固定 Q缰儿，優(yōu)化 θ；交替將極值推向極大勺美。Kmeans也是這樣更新的叭喜，固定中心點(diǎn)着倾，優(yōu)化Ein；優(yōu)化完成泣棋，更新中心點(diǎn)修赞。SMO也是，固定更新α1和α2荔泳。

④GMM的推導(dǎo)

可以直接把高斯混合模型代入EM框架里面沾鳄。
存在多個高斯分布混合生成了一堆數(shù)據(jù)X吞歼，取各個高斯分布的概率是糯俗，第i個高斯分布的均值是淘正，方差是惩淳，求法φ，μ含蓉，σ差油。
按照套路刃鳄，第一個E-step求出Q，于是有：
意思就是求出第i個樣本屬于第j個分布的概率是多少。之后就是M-step了伦忠，就是化簡了：

這里可能需要解釋一下，根據(jù)至于條件，因?yàn)楹苊黠@赋咽，z是隱變量旧噪，只是指明了x是屬于哪個類別，和μ脓匿，Σ沒有什么關(guān)系淘钟，所以直接忽略那兩個參數(shù)了，所以P(z)是沒有那兩個參數(shù)的陪毡，z是代表了分布米母，所以每一個分布的概率肯定是包括了，所以就只有一個概率的參數(shù)毡琉。P(x|z)是本身的概率铁瞒，就是已經(jīng)知道分布是那個了，求屬于這個分布的概率是多少桅滋，既然已經(jīng)選定了分布那么自然就不需要再看φ了慧耍，因?yàn)棣帐歉鱾€分布的概率。

⑤硬幣模型

現(xiàn)在有兩個硬幣AB擎勘，進(jìn)行5次試驗(yàn)每一次投10次，并不知道是哪個硬幣投的颖榜，求兩種硬幣的正面的概率棚饵。
首先E-step：
首先先初始化一下煤裙，
第一個試驗(yàn)選中A的概率：
同樣求得
計(jì)算機(jī)出每一個試驗(yàn)的概率然后相加求均值。
之后就是M-step了：

⑥代碼實(shí)現(xiàn)

方差的求解就不玩了欣硼，主要就是迭代求解μ和φ的值了题翰。
首先是生成數(shù)據(jù)，4個高斯分布分别，每一個高斯分布的sigma都是一樣的遍愿，不一樣的只有μ和α，也就是φ耘斩，習(xí)慣上把前面的一個參數(shù)叫做權(quán)值沼填，所以用α來表示。

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來找到他們屬于的類別坞笙。

其實(shí)如果用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即可薛夜，求解每一個分布對樣本點(diǎn)做了多少的功，之后求單個樣本點(diǎn)求比例版述。
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

直接按照公式優(yōu)化即可梯澜。

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)

運(yùn)行函數(shù)】饰觯看看結(jié)果：

結(jié)果其實(shí)還是相差不大晚伙。達(dá)到預(yù)期。

⑦EM的另一種理解方法

上面所講的其實(shí)只是一種理解方法俭茧，在李航老師的統(tǒng)計(jì)學(xué)習(xí)方法里面是另一種比較厲害的解法：
1.E-step：求出Q函數(shù)咆疗。
2.M-step：利用Q函數(shù)求極大值。
其實(shí)這兩種方法是完全一樣的母债，Q函數(shù)就是下界函數(shù)午磁，

這個公式和上面的：

是一樣的，至于為什么有下面這個

這是李航老師書上的毡们，既然θ是已經(jīng)知道了迅皇，那么自然就可以去掉，因?yàn)樵趌og下面做分母就相當(dāng)于是常數(shù)了衙熔，拆開可以變減號嘛喧半，在前面的就不可以∏嘣穑回來看看我們剛剛推導(dǎo)的挺据，其實(shí)是一樣的，下面的那個Q確實(shí)可以去掉脖隶，因?yàn)槭鞘孪戎赖牧吮饽汀Ｔ谑褂肑ensen不等式的時候，需要假設(shè)隱變量服從某種形式的概率分布产阱，才可以將推導(dǎo)過程的一部分看成是期望的表達(dá)形式從而應(yīng)用Jensen不等式婉称。然而這個分布不是隨便指定的。我們令Jensen不等式取等號的時候构蹬，可以計(jì)算出這個分布其實(shí)就是：已知觀測數(shù)據(jù)的隱變量的后驗(yàn)概率分布王暗。由于求Q函數(shù)需要先求出隱變量的后驗(yàn)概率的期望，因此庄敛，這就可以解釋為什么EM算法的“通俗”理解角度的E步驟是求隱變量的期望了俗壹。有時候在用EM算法解決某個具體問題的時候，會發(fā)現(xiàn)M步驟極大化的居然是完全數(shù)據(jù)的對數(shù)似然函數(shù)藻烤。這是因?yàn)楸脸琎函數(shù)雖然是完全數(shù)據(jù)的對數(shù)似然函數(shù)的某種期望，但是求這個期望的過程有時其實(shí)就是將隱變量的后驗(yàn)概率的期望代入就可以了怖亭。因此涎显，本質(zhì)上我們其實(shí)還是在求Q函數(shù)的極大。

⑧總結(jié)

EM和Kmeans算法其實(shí)很類似，事實(shí)上步驟基本可以用EM框架來替換拯刁，但是Kmeans算法是硬分類脊岳，說一不二，但是EM算法不太一樣垛玻，是軟分類割捅，百分之幾是那個，百分之幾是這個帚桩。

缺點(diǎn)也還是有的：初值敏感亿驾，局部最優(yōu)。因?yàn)榇嬖诹穗[變量账嚎，所以導(dǎo)致了直接對x做極大似然是不可行的莫瞬，log已經(jīng)在sum的外面了儡蔓。所以EM算法就轉(zhuǎn)向了下界函數(shù)，而這種方法本來就不保證找到局部最優(yōu)解疼邀。

如果將樣本看作觀察值喂江，潛在類別看作是隱藏變量，那么聚類問題也就是參數(shù)估計(jì)問題旁振。如果一個目標(biāo)函數(shù)存在多個變量获询，那么梯度下降牛頓法這些逼近方法就用不了了。但我們可以使用坐標(biāo)上升方法拐袜，固定一個變量吉嚣，對另外一個求導(dǎo)數(shù)，然后替換最后逐步逼近極值點(diǎn)蹬铺。對應(yīng)到EM算法也是一樣尝哆，E步求隱含的z變量，Mstep求解其他參數(shù)丛塌。

最后符上GitHub代碼：

https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/ExpectationMaximization