前言：前一篇文章大概說了EM算法的整個理解以及一些相關(guān)的公式神馬的，那些數(shù)學公式啥的看完真的是忘完了，那就來用代碼記憶記憶吧！接下來將會對python版本的EM算法進行一些分析。

EM的python實現(xiàn)和解析

引入問題（雙硬幣問題）

假設(shè)有兩枚硬幣A、B，以相同的概率隨機選擇一個硬幣，進行如下的拋硬幣實驗：共做5次實驗，每次實驗獨立的拋十次，結(jié)果如圖中a所示，例如某次實驗產(chǎn)生了H、T、T、T、H、H、T、H、T、H，H代表正面朝上。

假設(shè)試驗數(shù)據(jù)記錄員可能是實習生，業(yè)務(wù)不一定熟悉，造成a和b兩種情況

a表示實習生記錄了詳細的試驗數(shù)據(jù)，我們可以觀測到試驗數(shù)據(jù)中每次選擇的是A還是B

b表示實習生忘了記錄每次試驗選擇的是A還是B，我們無法觀測實驗數(shù)據(jù)中選擇的硬幣是哪個

問在兩種情況下分別如何估計兩個硬幣正面出現(xiàn)的概率？

以上的針對于b實習生的問題其實和三硬幣問題類似，只是這里把三硬幣中第一個拋硬幣的選擇換成了實習生的選擇。

對于已知是A硬幣還是B硬幣拋出的結(jié)果的時候，可以直接采用概率的求法來進行求解。對于含有隱變量的情況，也就是不知道到底是A硬幣拋出的結(jié)果還是B硬幣拋出的結(jié)果的時候，就需要采用EM算法進行求解了。如下圖：

其中的EM算法的第一步就是初始化的過程，然后根據(jù)這個參數(shù)得出應(yīng)該產(chǎn)生的結(jié)果。

構(gòu)建觀測數(shù)據(jù)集

針對這個問題，首先采集數(shù)據(jù)，用1表示H（正面），0表示T（反面）：

#硬幣投擲結(jié)果observations = numpy.array([[1,0,0,0,1,1,0,1,0,1],            [1,1,1,1,0,1,1,1,0,1],            [1,0,1,1,1,1,1,0,1,1],            [1,0,1,0,0,0,1,1,0,0],            [0,1,1,1,0,1,1,1,0,1]])

第一步：參數(shù)的初始化

參數(shù)賦初值

第一個迭代的E步

拋硬幣是一個二項分布，可以用scipy中的binom來計算。對于第一行數(shù)據(jù)，正反面各有5次，所以：

#二項分布求解公式contribution_A = scipy.stats.binom.pmf(num_heads,len_observation,theta_A)contribution_B = scipy.stats.binom.pmf(num_heads,len_observation,theta_B)

將兩個概率正規(guī)化，得到數(shù)據(jù)來自硬幣A，B的概率：

weight_A = contribution_A / (contribution_A + contribution_B)weight_B = contribution_B / (contribution_A + contribution_B)

這個值類似于三硬幣模型中的μ，只不過多了一個下標，代表是第幾行數(shù)據(jù)（數(shù)據(jù)集由5行構(gòu)成）。同理，可以算出剩下的4行數(shù)據(jù)的μ。