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

EM的python實現和解析

引入問題(雙硬幣問題)

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

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

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

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

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

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

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

其中的EM算法的第一步就是初始化的過程，然后根據這個參數得出應該產生的結果。

構建觀測數據集

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

#硬幣投擲結果

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]])

第一步：參數的初始化

參數賦初值

第一個迭代的E步

拋硬幣是一個二項分布，可以用scipy中的binom來計算。對于第一行數據，正反面各有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)

將兩個概率正規化，得到數據來自硬幣A，B的概率：

weight_A = contribution_A / (contribution_A + contribution_B)

weight_B = contribution_B / (contribution_A + contribution_B)

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

有了μ，就可以估計數據中AB分別產生正反面的次數了。μ代表數據來自硬幣A的概率的估計，將它乘上正面的總數，得到正面來自硬幣A的總數，同理有反面，同理有B的正反面。

#更新在當前參數下A，B硬幣產生的正反面次數

counts['A']['H'] += weight_A * num_heads

counts['A']['T'] += weight_A * num_tails

counts['B']['H'] += weight_B * num_heads

counts['B']['T'] += weight_B * num_tails

第一個迭代的M步

當前模型參數下，AB分別產生正反面的次數估計出來了，就可以計算新的模型參數了：

new_theta_A = counts['A']['H']/(counts['A']['H'] + counts['A']['T'])

new_theta_B = counts['B']['H']/(counts['B']['H'] + counts['B']['T'])

于是就可以整理一下，給出EM算法單個迭代的代碼：

def em_single(priors,observations):

"""

EM算法的單次迭代

Arguments

------------

priors:[theta_A,theta_B]

observation:[m X n matrix]

Returns

---------------

new_priors:[new_theta_A,new_theta_B]

:param priors:

:param observations:

:return:

"""

counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}}

theta_A = priors[0]

theta_B = priors[1]

#E step

for observation in observations:

len_observation = len(observation)

num_heads = observation.sum()

num_tails = len_observation-num_heads

#二項分布求解公式

contribution_A = scipy.stats.binom.pmf(num_heads,len_observation,theta_A)

contribution_B = scipy.stats.binom.pmf(num_heads,len_observation,theta_B)

weight_A = contribution_A / (contribution_A + contribution_B)

weight_B = contribution_B / (contribution_A + contribution_B)

#更新在當前參數下A，B硬幣產生的正反面次數

counts['A']['H'] += weight_A * num_heads

counts['A']['T'] += weight_A * num_tails

counts['B']['H'] += weight_B * num_heads

counts['B']['T'] += weight_B * num_tails

# M step

new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])

new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])

return [new_theta_A,new_theta_B]

EM算法主循環

給定循環的兩個終止條件：模型參數變化小于閾值；循環達到最大次數，就可以寫出EM算法的主循環了

def em(observations,prior,tol = 1e-6,iterations=10000):

"""

EM算法

：param observations :觀測數據

：param prior：模型初值

：param tol：迭代結束閾值

：param iterations：最大迭代次數

：return：局部最優的模型參數

"""

iteration = 0;

while iteration < iterations:

new_prior = em_single(prior,observations)

delta_change = numpy.abs(prior[0]-new_prior[0])

if delta_change < tol:

break

else:

prior = new_prior

iteration +=1

return [new_prior,iteration]

調用

給定數據集和初值，就可以調用EM算法了：

print em(observations,[0.6,0.5])

得到

[[0.72225028549925996, 0.55543808993848298], 36]

我們可以改變初值，試驗初值對EM算法的影響。

print em(observations,[0.5,0.6])

結果：

[[0.55543727869042425, 0.72225099139214621], 37]

看來EM算法還是很健壯的。如果把初值設為相等會怎樣？

print em(observations,[0.3,0.3])

輸出：[[0.64000000000000001, 0.64000000000000001], 1]

顯然，兩個值相加不為1的時候就會破壞這個EM函數。

換一下初值：

print em(observations,[0.99999,0.00001])

輸出：[[0.72225606292866507, 0.55543145006184214], 33]

EM算法對于參數的改變還是有一定的健壯性的。

以上是根據前人寫的博客進行學習的~可以自己動手實現以下，對于python練習還是有作用的。希望對大家的學習有所幫助，也希望大家多多支持。

總結

以上是生活随笔為你收集整理的python实现em聚类算法_EM算法的python实现的方法步骤的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。