現在有兩個硬幣a和b,要估計的引數是它們各自翻正面(head)的概率。觀察的過程是先隨機選a或者b,然後扔10次。以上步驟重複5次。如果知道每次選的是a還是b,那可以直接估計(見下圖a)。如果不知道選的是a還是b(隱變數),只觀測到5次迴圈共50次投幣的結果,這時就沒法直接估計a和b的正面概率。em演算法此時可起作用(見下圖b)。
稍微解釋一下上圖的計算過程。初始值θa=0.6,θb=0.5。
圖中的0.45是怎麼得來的呢?由兩個硬幣的初始值0.6和0.5,容易得出投擲出5正5反的概率是pa=c(10,5)(0.65)*(0.45),pb=c(10,5)(0.55)*(0.55), pa/(pa+pb)=0.449, 0.45就是0.449近似而來的,表示第一組實驗選擇的硬幣是a的概率為0.45。圖中的2.2h,2.2t是怎麼得來的呢? 0.449 * 5h = 2.2h ,0.449 * 5t = 2.2t ,表示第一組實驗選擇a硬幣且正面朝上次數的期望值是2.2。其他的值依次類推。
python**實現:
import numpy as np
import math
input = np.array([[1, 0, 0, 0, 1, 1, 0, 1, 0, 1],
[1, 1, 1, 1, 0, 1, 1, 1, 1, 1],
[1, 1, 0, 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]])
thetaa = 0.5
thetab = 0.6
m, n = input.shape
headnums = input.sum(axis = 1)
tailnums = n - headnums
for i in range(10):
pahnum = 0
patnum = 0
pbhnum = 0
pbtnum = 0
for headnum in headnums:
pa = math.pow(thetaa, headnum) * math.pow((1 - thetaa), (n - headnum))
pb = math.pow(thetab, headnum) * math.pow((1 - thetab), (n - headnum))
psum = pa + pb
pa = pa / psum
pb = pb / psum
pahnum += headnum * pa
patnum += (n - headnum)*pa
pbhnum += headnum * pb
pbtnum += (n - headnum)*pb
thetaa = pahnum / (pahnum + patnum)
thetab = pbhnum / (pbhnum + pbtnum)
EM 演算法 例項
coding utf 8 import math import copy import numpy as np import matplotlib.pyplot as plt isdebug true 指定k個高斯分布引數,這裡指定k 2。注意2個高斯分布具有相同方差sigma,均值分別為mu1,m...
簡單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演算法之Python
之 示例一 二硬幣模型 假設現在有兩個硬幣a和b,我們想要知道兩枚硬幣各自為正面的概率啊即模型的引數。我們先隨機從a,b中選一枚硬幣,然後扔10次並記錄下相應的結果,h代表正面t代表反面。對以上的步驟重複進行5次。如果在記錄的過程中我們記錄下來每次是哪一枚硬幣 即知道每次選的是a還是b 那可以直接根...