您好,登錄后才能下訂單哦!
今天就跟大家聊聊有關(guān)使用Python實(shí)現(xiàn)EM算法,可能很多人都不太了解,為了讓大家更加了解,小編給大家總結(jié)了以下內(nèi)容,希望大家根據(jù)這篇文章可以有所收獲。
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ù)θ就非常簡單了,如
下圖所示:
如果不知道每次拋的是哪個(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))
看完上述內(nèi)容,你們對(duì)使用Python實(shí)現(xiàn)EM算法有進(jìn)一步的了解嗎?如果還想了解更多知識(shí)或者相關(guān)內(nèi)容,請(qǐng)關(guān)注億速云行業(yè)資訊頻道,感謝大家的支持。
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng),如果涉及侵權(quán)請(qǐng)聯(lián)系站長郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。