**之：
示例一：二硬幣模型
假設現在有兩個硬幣a和b，我們想要知道兩枚硬幣各自為正面的概率啊即模型的引數。我們先隨機從a,b中選一枚硬幣，然後扔10次並記錄下相應的結果，h代表正面t代表反面。對以上的步驟重複進行5次。如果在記錄的過程中我們記錄下來每次是哪一枚硬幣（即知道每次選的是a還是b），那可以直接根據結果進行估計（見下圖a）。
但是如果資料中沒記錄每次投擲的硬幣是a還是b（隱變數），只觀測到5次迴圈共50次投幣的結果，這時就沒法直接估計a和b的正面概率。這時就該輪到em演算法大顯身手了，em演算法特別適用於這種含有隱變數的引數求解問題（見下圖b）。
示例二：三硬幣模型
現在我們將上面的二硬幣模型擴充套件為三硬幣模型，其實原理基本差不多。假設有三枚硬幣a、b、c，這些硬幣正面出現的概率分別p，q和
。先拋c硬幣，如果c硬幣為正面則選擇硬幣a，反之選擇硬幣b，然後對選出的硬幣進行一組實驗，獨立的拋十次。共做5次實驗，每次實驗獨立的拋十次，結果如圖中a所示，例如某次實驗產生了h、t、t、t、h、h、t、h、t、h，h代表正面朝上。
本人最近也剛學em演算法，下面**主要參考em演算法及其推廣，這裡面作者實現了乙個兩硬幣模型的em演算法。本文對其稍做了一點修改，變成三硬幣模型。
em演算法步驟：
e步：計算在當前迭代的模型引數下，觀測資料y來自硬幣b的概率：
對於這個三硬幣模型來說，我們先通過e步（對隱變數求期望）來求得隱變數的引數（即屬於哪個硬幣），然後再通過m-step來重新估算三個硬幣的引數，直至收斂（達到要求）為止。下面是實現三硬幣模型的em演算法**，希望可以更好的幫助理解。
#執行後結果為：!usr/bin/env python
#-*- coding:utf-8 -*-
import
numpy as np
from scipy import
stats
defem_single(priors, observations):
"""em演算法單次迭代
arguments
---------
priors : [theta_a, theta_b，theta_c]
observations : [m x n matrix]
returns
--------
new_priors: [new_theta_a, new_theta_b,new_theta_c]
:param priors:
:param observations:
:return:
"""counts = , '
b': }
theta_a =priors[0]
theta_b = priors[1]
theta_c=priors[2]
#e step
weight_as=
for observation in
observations:
len_observation =len(observation)
num_heads =observation.sum()
num_tails = len_observation -num_heads
contribution_a = theta_c*stats.binom.pmf(num_heads, len_observation, theta_a)
contribution_b = (1-theta_c)*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_c = 1.0*sum(weight_as)/len(weight_as)
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,new_theta_c]
def em(observations, prior, tol=1e-6, iterations=10000):
"""em演算法
:param observations: 觀測資料
:param prior: 模型初值
:param tol: 迭代結束閾值
:param iterations: 最大迭代次數
:return: 區域性最優的模型引數
"""import
math
iteration =0
while iteration
new_prior =em_single(prior, observations)
delta_change = np.abs(prior[0] -new_prior[0])
if delta_change
break
else
: prior =new_prior
iteration += 1
return
[new_prior, iteration]
#硬幣投擲結果觀測序列：1表示正面，0表示反面。
observations = np.array([[1, 0, 0, 0, 1, 1, 0, 1, 0, 1],
[1, 1, 1, 1, 0, 1, 1, 1, 1, 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.5, 0.8, 0.6])
[[0.51392121603987106, 0.79337052912023864, 0.47726196801164544], 42]
從結果我們可以了解到經過42輪迭代，我們最終得出了結果：硬幣a正面的概率為0.51392121603987106，硬幣b為正面的概率為0.79337052912023864，c硬幣正面概率為0.47726196801164544。
至此em演算法的實現就完成了
簡單EM演算法python實現
reference coding utf 8 from future import print function 知道投擲的是a或b,求coin a和b投擲結果是head的最大似然估計,用來評估em演算法 mle a 24.0 24 6 mle b 9.0 9 11 print real mle a...
燈管實驗的em演算法 EM演算法
本文試圖用最簡單的例子 最淺顯的方式說明em expectation maximization 演算法的應用場景和使用方法，而略去公式的推導和收斂性的證明。maximum likelihood estimation maximum likelihood estimation mle 是要選擇乙個最佳...
EM演算法 存在雜訊的EM演算法
1 em演算法的 r語言實現 步驟1 資料集準備及其描述 步驟2 構建em 演算法模型，指定分3類 步驟3 構建em 演算法模型，指定先驗概率 2 存在雜訊的 em演算法r語言實現 步驟1 資料集準備 set.seed 0 設定隨機種子 nnoise 100 dim faithful 解釋 runi...