溫馨提示×
您好，登錄后才能下訂單哦！
點擊 登錄注冊 即表示同意《億速云用戶服務條款》
您好，登錄后才能下訂單哦！
今天就跟大家聊聊有關使用Python實現EM算法，可能很多人都不太了解，為了讓大家更加了解，小編給大家總結了以下內容，希望大家根據這篇文章可以有所收獲。
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實現
#coding=utf-8
from numpy import *
from scipy import 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 observation in 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算法有進一步的了解嗎？如果還想了解更多知識或者相關內容，請關注億速云行業資訊頻道，感謝大家的支持。
免責聲明：本站發布的內容（圖片、視頻和文字）以原創、轉載和分享為主，文章觀點不代表本網站立場，如果涉及侵權請聯系站長郵箱：is@yisu.com進行舉報，并提供相關證據，一經查實，將立刻刪除涉嫌侵權內容。
