參考：

一、一個簡單的概率問題

實驗：現在有A和B兩個硬幣，我們從這兩個硬幣中，隨機選取5次，做5組實驗，每組實驗內容是：丟擲所選的硬幣，記錄正反面。

實驗資料如下：

目標：根據所得到的實驗資料，分別求出硬幣A和B丟擲後正面向上的概率。

根據古典概率的原理，容易求出:

硬幣A正面向上的概率為 PA = (9 + 8 + 7) / (10*3) = 0.8

硬幣B正面向上的概率為 PB = (5 + 4) / (10*2) = 0.45

二、如果不知道所選擇的的硬幣呢

如果不知道所選取的硬幣是A還是B，只記錄硬幣拋下後的正反面，那麼記錄的資料如下：

在這種情況下，目標相同：根據上面的實驗資料，分別求出硬幣A和B丟擲後正面向上的概率。

EM演算法的做法如下：

隨機設定PA = 0.8, PB = 0.6，我們計算一下，針對每次實驗，我們分別假設是硬幣A還是B，然後計算出現對應實驗結果的概率。比如針對第一次實驗，假設選擇的硬幣是A，那麼出現對應的實驗結果（5正5反）的概率為：

0.8^5*(1-0.8)^5

如此，計算對應的概率分別為：

反推的結果意思是，比較每組的兩個概率。比如，針對第一次實驗，選擇A得到的實驗結果的概率<B，所以反推得到本次所選擇的硬幣為B。

因此，計算PA和PB:

PA = (9 + 8) / (10 * 2) = 0.85

PB = (5 + 4 + 7) / (10 * 3) = 0.5333...

可以看出這和之前假設的PA = 0.85， PB = 0.6有所不同。這個過程是一個迭代的過程，一直用上面的方法，更新PA和PB，直到收斂，得到的PA和PB就是最終的答案。

上述問題的Python程式碼如下：

import numpy as np


if __name__ == '__main__':

    n, m = 5, 10    # 5組實驗，每組10次
    nzs = [5, 9, 8, 4, 7]   # 每組實驗，硬幣正面的次數

    # 初始化
    PA = 0.8    # 待求的概率：拋硬幣A正面向上的概率
    PB = 0.6    # 待求的概率：拋硬幣B正面向上的概率

    n_iter = 10
    for i in range(0, n_iter):  # 迭代多次
        selected_coins = []
        for j in range(0, n):
            # P_A 記錄的是，如果第j次實驗，選擇的是硬幣A，那麼出現對應的實驗結果的概率
            P_A = np.power(PA, nzs[j])*np.power(1-PA, m-nzs[j])
            P_B = np.power(PB, nzs[j])*np.power(1-PB, m-nzs[j])
            if P_A > P_B:
                selected_coins.append('A')
            else:
                selected_coins.append('B')
        fenzi = [0, 0]
        fenmu = [0, 0]
        for j, x in enumerate(selected_coins):
            if x == 'A':
                fenzi[0] += nzs[j]
                fenmu[0] += m
            else:
                fenzi[1] += nzs[j]
                fenmu[1] += m

        PA = fenzi[0] / fenmu[0]
        PB = fenzi[1] / fenmu[1]
        print(PA, PB)
 

迭代10次之後，執行結果如下：

可以看出，第二次就收斂了。真實的值是PA = 0.8, PB = 0.45。EM演算法求得結果與實際結果有差別。

三、小結

上述過程，可以歸納為如下幾個步驟：

step1. 隨機設定PA和PB的值

step2. 反推所取的硬幣是A還是B（Maximum）

step3. 根據反推的結果重新計算PA和PB（Expectation）

step4. 迭代上述過程

以上就是期望最大化的過程。其中step2，“反推”的過程採用的是極大似然估計。

四、改進的方法

上述案例中，step2“反推”的過程採用極大似然估計，“反推”所選取的硬幣是A還是B，改進的方法是：將概率歸一化如下：

針對第一組資料，此時不再根據 0.116 < 0.884，直接判定，第一組選擇的硬幣是B，而是說，選擇硬幣A的概率是0.116，選擇硬幣B的概率是0.884。

在這種情況下，PA的計算方法如下：

PA = 正面的情況 / [正面的情況 + (反面的情況)] 

正面的情況 = (5*0.116 + 9*0.869 + 8*0.714 + 4*0.0471 + 7*0.484)

反面的情況 = (5*0.116 + 1*0.869 + 2*0.714 + 4(0.0471 + 3*0.484)

正面的情況 + 反面的情況 = 10 * （0.116 + 0.869 + 0.714 + 0.0471 + 0.484）

對應的程式碼如下：

import numpy as np


if __name__ == '__main__':

    n, m = 5, 10    # 5組實驗，每組10次
    nzs = [5, 9, 8, 4, 7]   # 每組實驗，硬幣正面的次數

    # 初始化
    PA = 0.8    # 待求的概率：拋硬幣A正面向上的概率
    PB = 0.6    # 待求的概率：拋硬幣B正面向上的概率

    n_iter = 10
    for i in range(0, n_iter):  # 迭代多次
        selected_coins_p_a = []
        selected_coins_p_b = []
        for j in range(0, n):
            # P_A 記錄的是，如果第j次實驗，選擇的是硬幣A，那麼出現對應的實驗結果的概率
            P_A = np.power(PA, nzs[j])*np.power(1-PA, m-nzs[j])
            P_B = np.power(PB, nzs[j])*np.power(1-PB, m-nzs[j])

            # 歸一化
            P_A = P_A / (P_A + P_B)
            P_B = 1 - P_A
            
            selected_coins_p_a.append(P_A)
            selected_coins_p_b.append(P_B)
            #
            # print('P_A', P_A)
            # print('P_B', P_B)
            # exit()

        zhengmian, fanmian = 0, 0
        for j in range(0, n):
            zhengmian += nzs[j] * selected_coins_p_a[j]
            fanmian += (m - nzs[j]) * selected_coins_p_a[j]
        PA = zhengmian / (zhengmian + fanmian)

        zhengmian, fanmian = 0, 0
        for j in range(0, n):
            zhengmian += nzs[j] * selected_coins_p_b[j]
            fanmian += (m - nzs[j]) * selected_coins_p_b[j]
        PB = zhengmian / (zhengmian + fanmian)

        print(PA, PB)

執行結果如下：

真實結果是0.8, 0.45。同樣是迭代10次，改進之前的結果是0.85, 0.53。改進之後的結果是0.79， 0.51，更接近真實值。

本文通過簡單的案例，理解了EM演算法背後的原理。