今天就跟大家聊聊有關(guān)使用Python實(shí)現(xiàn)EM算法，可能很多人都不太了解，為了讓大家更加了解，小編給大家總結(jié)了以下內(nèi)容，希望大家根據(jù)這篇文章可以有所收獲。

自2013年創(chuàng)立以來創(chuàng)新互聯(lián)建站專注于”幫助中小企業(yè)+互聯(lián)網(wǎng)”, 也是目前成都地區(qū)具有實(shí)力的互聯(lián)網(wǎng)服務(wù)商。團(tuán)隊(duì)致力于為企業(yè)提供--站式網(wǎng)站建設(shè)、移動(dòng)端應(yīng)用( H5手機(jī)營銷、重慶APP開發(fā)公司、微信開發(fā))、軟件開發(fā)、信息化解決方案等服務(wù)。

EM算法實(shí)例

通過實(shí)例可以快速了解EM算法的基本思想，具體推導(dǎo)請(qǐng)點(diǎn)文末鏈接。圖a是讓我們預(yù)熱的，圖b是EM算法的實(shí)例。

這是一個(gè)拋硬幣的例子，H表示正面向上，T表示反面向上，參數(shù)θ表示正面朝上的概率。硬幣有兩個(gè)，A和B，硬幣是有偏的。本次實(shí)驗(yàn)總共做了5組，每組隨機(jī)選一個(gè)硬幣，連續(xù)拋10次。如果知道每次拋的是哪個(gè)硬幣，那么計(jì)算參數(shù)θ就非常簡(jiǎn)單了，如

下圖所示：

如果不知道每次拋的是哪個(gè)硬幣呢？那么，我們就需要用EM算法，基本步驟為：

  1、給θ_AθA​和θ_BθB​一個(gè)初始值；

  2、（E-step）估計(jì)每組實(shí)驗(yàn)是硬幣A的概率（本組實(shí)驗(yàn)是硬幣B的概率=1-本組實(shí)驗(yàn)是硬幣A的概率）。分別計(jì)算每組實(shí)驗(yàn)中，選擇A硬幣且正面朝上次數(shù)的期望值，選擇B硬幣且正面朝上次數(shù)的期望值；

  3、（M-step）利用第三步求得的期望值重新計(jì)算θ_AθA​和θ_BθB​；

  4、當(dāng)?shù)揭欢ù螖?shù)，或者算法收斂到一定精度，結(jié)束算法，否則，回到第2步。

計(jì)算過程詳解：初始值θ_A^{(0)}θA(0)​=0.6,θ_B^{(0)}θB(0)​=0.5。

由兩個(gè)硬幣的初始值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近似而來的，表示第一組實(shí)驗(yàn)選擇的硬幣是A的概率為0.45。然后，0.449 * 5H = 2.2H ，0.449 * 5T = 2.2T ，表示第一組實(shí)驗(yàn)選擇A硬幣且正面朝上次數(shù)和反面朝上次數(shù)的期望值都是2.2，其他的值依次類推。最后，求出θ_A^{(1)}θA(1)​=0.71,θ_B^{(1)}θB(1)​=0.58。重復(fù)上述過程，不斷迭代，直到算法收斂到一定精度為止。

這篇博客對(duì)EM算法的推導(dǎo)非常詳細(xì)，鏈接如下：

https://blog.csdn.net/zhihua_oba/article/details/73776553

Python實(shí)現(xiàn)

#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
  #二項(xiàng)分布求解公式
  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)
  #更新在當(dāng)前參數(shù)下A，B硬幣產(chǎn)生的正反面次數(shù)
  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 :觀測(cè)數(shù)據(jù)
 ：param prior：模型初值
 ：param tol：迭代結(jié)束閾值
 ：param iterations：大迭代次數(shù)
 ：return：局部最優(yōu)的模型參數(shù)
 """
 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]

#硬幣投擲結(jié)果
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實(shí)現(xiàn)EM算法-創(chuàng)新互聯(lián)
標(biāo)題路徑：http://www.rwnh.cn/article8/epiip.html

成都網(wǎng)站建設(shè)公司_創(chuàng)新互聯(lián)，為您提供網(wǎng)站收錄、動(dòng)態(tài)網(wǎng)站、App開發(fā)、微信小程序、網(wǎng)站設(shè)計(jì)、網(wǎng)站制作

廣告

聲明：本網(wǎng)站發(fā)布的內(nèi)容（圖片、視頻和文字）以用戶投稿、用戶轉(zhuǎn)載內(nèi)容為主，如果涉及侵權(quán)請(qǐng)盡快告知，我們將會(huì)在第一時(shí)間刪除。文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng)，如需處理請(qǐng)聯(lián)系客服。電話：028-86922220；郵箱：631063699@qq.com。內(nèi)容未經(jīng)允許不得轉(zhuǎn)載，或轉(zhuǎn)載時(shí)需注明來源： 創(chuàng)新互聯(lián)