#coding:utf-8
import math
import copy
import numpy as np
import matplotlib.pyplot as plt
isdebug = true
#指定k個高斯分布引數，這裡指定k=2。
#注意2個高斯分布具有相同方差sigma，均值分別為mu1,mu2。
#共1000個資料
#生成訓練樣本，輸入6,40,20,2
#兩類樣本方差為6，
#一類均值為20，一類為40
#隨機生成1000個數
def ini_data(sigma,mu1,mu2,k,n):
#儲存生成的隨機樣本
global x
#求類別的均值
global mu
#儲存樣本屬於某類的概率
global expectations
#1*n的矩陣，生成n個樣本
x = np.zeros((1,n))
#任意給定兩個初始值，任猜兩類均值
#賦值一次即可，最後要輸出的量
mu = np.random.random(2) #0-1之間
print mu
#給定1000*2的矩陣，儲存樣本屬於某類的概率
expectations = np.zeros((n,k))
#生成n個樣本資料
for i in xrange(0,n):
#在大於0.5在第1個分布，小於0.5在第2個分布
if np.random.random(1) > 0.5:
#均值40加上方差倍數，樣本資料滿足n(40,sigma)正態分佈
x[0,i] = np.random.normal()*sigma + mu1 #
else:
#均值40加上方差倍數，樣本資料滿足n(20,sigma)正態分佈
x[0,i] = np.random.normal()*sigma + mu2
if isdebug:
print "***********"
print u"初始觀測資料x："
print x
#e步 計算每個樣本屬於男女各自的概率
#輸入：方差sigma，類別k，樣本數n
def e_step(sigma,k,n):
#樣本屬於某類概率
global expectations
#兩類均值
global mu
#樣本global x
#遍歷所有樣本點，計算屬於每個類別的概率
for i in xrange(0,n):
#分母，用於歸一化
denom = 0
#遍歷男女兩類，計算各自歸一化分母
for j in xrange(0,k):
#計算分母
denom += math.exp((-1/(2*(float(sigma**2))))*(float(x[0,i]-mu[j]))**2)
#遍歷男女兩類，計算各自分子部分
for j in xrange(0,k):
#分子numer = math.exp((-1/(2*(float(sigma**2))))*(float(x[0,i]-mu[j]))**2)
#每個樣本屬於該類別的概率
expectations[i,j] = numer/denom
if isdebug:
print "***********"
print u"隱藏變數e（z）："
print len(expectations)
#資料總個數
print expectations.size
#矩陣資料
print expectations.shape
#列印出隱藏變數的值
print expectations
#m步 期望最大化
def m_step(k,n):
#樣本屬於某類概率p(k|xi)
global expectations
#樣本global x
#計算兩類的均值
#遍歷兩類
for j in xrange(0,k):
numer = 0
denom = 0
#當前類別下，遍歷所有樣本
#計算該類別下的均值和方差
for i in xrange(0,n):
#該類別樣本分佈p(k|xi)xi
numer += expectations[i,j]*x[0,i]
#該類別類樣本總數nk，nk等於p(k|xi)求和
denom +=expectations[i,j]
#計算每個類別各自均值uk
mu[j] = numer / denom
#演算法迭代iter_num次，或達到精度epsilon停止迭代
#迭代次數1000次， 誤差達到0.0001終止
#輸入：兩類相同方差sigma，一類均值mu1，一類均值mu2
#類別數k，樣本數n，迭代次數iter_num，可接受精度epsilon
def run(sigma,mu1,mu2,k,n,iter_num,epsilon):
#生成訓練樣本
ini_data(sigma,mu1,mu2,k,n)
print u"初始:", mu
#迭代1000次
for i in range(iter_num):
#儲存上次兩類均值
old_mu = copy.deepcopy(mu)
#e步e_step(sigma,k,n)
#m步m_step(k,n)
#輸出當前迭代次數及當前估計的值
print i,mu
#判斷誤差
if sum(abs(mu-old_mu)) < epsilon:
break
if __name__ == '__main__':
#sigma,mu1,mu2,模型數，樣本總數，迭代次數，迭代終止收斂精度
run(6,40,20,2,1000,1000,0.0001)
plt.hist(x[0,:],100) #柱狀圖的寬度
plt.show()
EM演算法例項及python實現
現在有兩個硬幣a和b，要估計的引數是它們各自翻正面 head 的概率。觀察的過程是先隨機選a或者b，然後扔10次。以上步驟重複5次。如果知道每次選的是a還是b，那可以直接估計 見下圖a 如果不知道選的是a還是b 隱變數 只觀測到5次迴圈共50次投幣的結果，這時就沒法直接估計a和b的正面概率。em演算...
燈管實驗的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...