EM算法的python實(shí)現(xiàn)的方法步驟
前言：前一篇文章大概說(shuō)了EM算法的整個(gè)理解以及一些相關(guān)的公式神馬的，那些數(shù)學(xué)公式啥的看完真的是忘完了，那就來(lái)用代碼記憶記憶吧！接下來(lái)將會(huì)對(duì)python版本的EM算法進(jìn)行一些分析。
EM的python實(shí)現(xiàn)和解析
引入問(wèn)題（雙硬幣問(wèn)題）
假設(shè)有兩枚硬幣A、B，以相同的概率隨機(jī)選擇一個(gè)硬幣，進(jìn)行如下的拋硬幣實(shí)驗(yàn)：共做5次實(shí)驗(yàn)，每次實(shí)驗(yàn)獨(dú)立的拋十次，結(jié)果如圖中a所示，例如某次實(shí)驗(yàn)產(chǎn)生了H、T、T、T、H、H、T、H、T、H，H代表正面朝上。
假設(shè)試驗(yàn)數(shù)據(jù)記錄員可能是實(shí)習(xí)生，業(yè)務(wù)不一定熟悉，造成a和b兩種情況
a表示實(shí)習(xí)生記錄了詳細(xì)的試驗(yàn)數(shù)據(jù)，我們可以觀測(cè)到試驗(yàn)數(shù)據(jù)中每次選擇的是A還是B
b表示實(shí)習(xí)生忘了記錄每次試驗(yàn)選擇的是A還是B，我們無(wú)法觀測(cè)實(shí)驗(yàn)數(shù)據(jù)中選擇的硬幣是哪個(gè)
問(wèn)在兩種情況下分別如何估計(jì)兩個(gè)硬幣正面出現(xiàn)的概率？
以上的針對(duì)于b實(shí)習(xí)生的問(wèn)題其實(shí)和三硬幣問(wèn)題類似，只是這里把三硬幣中第一個(gè)拋硬幣的選擇換成了實(shí)習(xí)生的選擇。
對(duì)于已知是A硬幣還是B硬幣拋出的結(jié)果的時(shí)候，可以直接采用概率的求法來(lái)進(jìn)行求解。對(duì)于含有隱變量的情況，也就是不知道到底是A硬幣拋出的結(jié)果還是B硬幣拋出的結(jié)果的時(shí)候，就需要采用EM算法進(jìn)行求解了。如下圖：
其中的EM算法的第一步就是初始化的過(guò)程，然后根據(jù)這個(gè)參數(shù)得出應(yīng)該產(chǎn)生的結(jié)果。
構(gòu)建觀測(cè)數(shù)據(jù)集
針對(duì)這個(gè)問(wèn)題，首先采集數(shù)據(jù)，用1表示H（正面），0表示T（反面）：
#硬幣投擲結(jié)果
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]])
第一步：參數(shù)的初始化
參數(shù)賦初值
第一個(gè)迭代的E步
拋硬幣是一個(gè)二項(xiàng)分布，可以用scipy中的binom來(lái)計(jì)算。對(duì)于第一行數(shù)據(jù)，正反面各有5次，所以：
#二項(xiàng)分布求解公式 contribution_A = scipy.stats.binom.pmf(num_heads,len_observation,theta_A) contribution_B = scipy.stats.binom.pmf(num_heads,len_observation,theta_B)
將兩個(gè)概率正規(guī)化，得到數(shù)據(jù)來(lái)自硬幣A，B的概率：
weight_A = contribution_A / (contribution_A + contribution_B) weight_B = contribution_B / (contribution_A + contribution_B)
這個(gè)值類似于三硬幣模型中的μ，只不過(guò)多了一個(gè)下標(biāo)，代表是第幾行數(shù)據(jù)（數(shù)據(jù)集由5行構(gòu)成）。同理，可以算出剩下的4行數(shù)據(jù)的μ。
有了μ，就可以估計(jì)數(shù)據(jù)中AB分別產(chǎn)生正反面的次數(shù)了。μ代表數(shù)據(jù)來(lái)自硬幣A的概率的估計(jì)，將它乘上正面的總數(shù)，得到正面來(lái)自硬幣A的總數(shù)，同理有反面，同理有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
第一個(gè)迭代的M步
當(dāng)前模型參數(shù)下，AB分別產(chǎn)生正反面的次數(shù)估計(jì)出來(lái)了，就可以計(jì)算新的模型參數(shù)了：
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算法單個(gè)迭代的代碼：
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 = 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)
#更新在當(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]
EM算法主循環(huán)
給定循環(huán)的兩個(gè)終止條件：模型參數(shù)變化小于閾值；循環(huán)達(dá)到最大次數(shù)，就可以寫出EM算法的主循環(huán)了
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 = numpy.abs(prior[0]-new_prior[0])
if delta_change < tol:
break
else:
prior = new_prior
iteration +=1
return [new_prior,iteration]
調(diào)用
給定數(shù)據(jù)集和初值，就可以調(diào)用EM算法了：
print em(observations,[0.6,0.5])
得到
[[0.72225028549925996, 0.55543808993848298], 36]
我們可以改變初值，試驗(yàn)初值對(duì)EM算法的影響。
print em(observations,[0.5,0.6])
結(jié)果：
[[0.55543727869042425, 0.72225099139214621], 37]
看來(lái)EM算法還是很健壯的。如果把初值設(shè)為相等會(huì)怎樣？
print em(observations,[0.3,0.3])
輸出：[[0.64000000000000001, 0.64000000000000001], 1]
顯然，兩個(gè)值相加不為1的時(shí)候就會(huì)破壞這個(gè)EM函數(shù)。
換一下初值：
print em(observations,[0.99999,0.00001])
輸出：[[0.72225606292866507, 0.55543145006184214], 33]
EM算法對(duì)于參數(shù)的改變還是有一定的健壯性的。
以上是根據(jù)前人寫的博客進(jìn)行學(xué)習(xí)的~可以自己動(dòng)手實(shí)現(xiàn)以下，對(duì)于python練習(xí)還是有作用的。希望對(duì)大家的學(xué)習(xí)有所幫助，也希望大家多多支持腳本之家。
相關(guān)文章
這篇文章主要為大家詳細(xì)介紹了如何使用Python開發(fā)一個(gè)在線編輯器,文中的示例代碼講解詳細(xì),具有一定的借鑒價(jià)值,有需要的小伙伴可以了解一下2025-02-02
Pandas中Concat與Append的實(shí)現(xiàn)與區(qū)別小結(jié)
本文主要介紹了Pandas中Concat與Append的實(shí)現(xiàn)與區(qū)別小結(jié),文中通過(guò)示例代碼介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2023-11-11
詳解Python如何循環(huán)遍歷Numpy中的Array
Numpy是Python中常見的數(shù)據(jù)處理庫(kù)，是數(shù)據(jù)科學(xué)中經(jīng)常使用的庫(kù)。在本文中，我們將學(xué)習(xí)如何迭代遍歷訪問(wèn)矩陣中的元素，需要的可以參考一下2022-04-04
Python+OpenCV手勢(shì)檢測(cè)與識(shí)別Mediapipe基礎(chǔ)篇
網(wǎng)上搜到了一些關(guān)于手勢(shì)處理的實(shí)驗(yàn),我在這兒簡(jiǎn)單的實(shí)現(xiàn)一下,下面這篇文章主要給大家介紹了關(guān)于Python+OpenCV手勢(shì)檢測(cè)與識(shí)別Mediapipe基礎(chǔ)篇的相關(guān)資料,需要的朋友可以參考下2022-12-12
Django-xadmin后臺(tái)導(dǎo)入json數(shù)據(jù)及后臺(tái)顯示信息圖標(biāo)和主題更改方式
這篇文章主要介紹了Django-xadmin后臺(tái)導(dǎo)入json數(shù)據(jù)及后臺(tái)顯示信息圖標(biāo)和主題更改方式，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2020-03-03
簡(jiǎn)單了解pytest測(cè)試框架setup和tearDown
這篇文章主要介紹了簡(jiǎn)單了解pytest測(cè)試框架setup和tearDown,文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2020-04-04
wxPython繪圖模塊wxPyPlot實(shí)現(xiàn)數(shù)據(jù)可視化
這篇文章主要為大家詳細(xì)介紹了wxPython繪圖模塊wxPyPlot實(shí)現(xiàn)數(shù)據(jù)可視化，文中示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2019-11-11
這篇文章主要為大家詳細(xì)介紹了python3爬取各類天氣信息，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2018-02-02
python忽略警告(warning)的3種方法小結(jié)
python開發(fā)中經(jīng)常遇到報(bào)錯(cuò)的情況,但是warning通常并不影響程序的運(yùn)行,而且有時(shí)特別討厭,下面我們來(lái)說(shuō)下如何忽略warning錯(cuò)誤,這篇文章主要給大家介紹了關(guān)于python忽略警告(warning)的3種方法,需要的朋友可以參考下2023-10-10
這篇文章主要介紹了python3 cmp實(shí)現(xiàn)方式，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。如有錯(cuò)誤或未考慮完全的地方，望不吝賜教2022-02-02