EM算法執行個體
通過執行個體可以快速了解EM算法的基本思想，具體推導請點文末連結。圖a是讓我們預熱的，圖b是EM算法的執行個體。
這是一個抛硬币的例子，H表示正面向上，T表示反面向上，參數θ表示正面朝上的機率。硬币有兩個，A和B，硬币是有偏的。本次實驗總共做了5組，每組随機選一個硬币，連續抛10次。如果知道每次抛的是哪個硬币，那麼計算參數θ就非常簡單了，如
下圖所示：
如果不知道每次抛的是哪個硬币呢？那麼，我們就需要用EM算法，基本步驟為：
1、給θ_AθA和θ_BθB一個初始值；
2、（E-step）估計每組實驗是硬币A的機率（本組實驗是硬币B的機率=1-本組實驗是硬币A的機率）。分别計算每組實驗中，選擇A硬币且正面朝上次數的期望值，選擇B硬币且正面朝上次數的期望值；
3、（M-step）利用第三步求得的期望值重新計算θ_AθA和θ_BθB；
4、當疊代到一定次數，或者算法收斂到一定精度，結束算法，否則，回到第2步。
計算過程詳解：初始值θ_A^{(0)}θA(0)=0.6,θ_B^{(0)}θB(0)=0.5。
由兩個硬币的初始值0.6和0.5，容易得出投擲出5正5反的機率是p_A=C^5_{10}*(0.6^5)*(0.4^5)pA=C105∗(0.65)∗(0.45)，p_B=C_{10}^5*(0.5^5)*(0.5^5)pB=C105∗(0.55)∗(0.55), p_ApA/(p_ApA+p_BpB)=0.449, 0.45就是0.449近似而來的，表示第一組實驗選擇的硬币是A的機率為0.45。然後，0.449 * 5H = 2.2H ，0.449 * 5T = 2.2T ，表示第一組實驗選擇A硬币且正面朝上次數和反面朝上次數的期望值都是2.2，其他的值依次類推。最後，求出θ_A^{(1)}θA(1)=0.71,θ_B^{(1)}θB(1)=0.58。重複上述過程，不斷疊代，直到算法收斂到一定精度為止。
這篇部落格對EM算法的推導非常詳細，連結如下：
https://blog.csdn.net/zhihua_oba/article/details/73776553
Python實作
?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#coding=utf-8
from numpyimport *
from scipyimport stats
import time
start= time.perf_counter()
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 observationin observations:
len_observation= len(observation)
num_heads= observation.sum()
num_tails= len_observation-num_heads
#二項分布求解公式
contribution_A= stats.binom.pmf(num_heads,len_observation,theta_A)
contribution_B= 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]
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= abs(prior[0]-new_prior[0])
if delta_change < tol:
break
else:
prior= new_prior
iteration+=1
return [new_prior,iteration]
#硬币投擲結果
observations= 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]])
print (em(observations,[0.6,0.5]))
end= time.perf_counter()
print('Running time: %f seconds'%(end-start))
總結
到此這篇關于Python實作EM算法執行個體的文章就介紹到這了,更多相關Python實作EM算法執行個體内容請搜尋伺服器之家以前的文章或繼續浏覽下面的相關文章希望大家以後多多支援伺服器之家！
原文連結：https://www.pianshen.com/article/7293314043/