#coding:gbk

import math

import copy

import numpy as np

import matplotlib.pyplot as plt

isdebug = False

# 指定k個高斯分布參數，這里指定k=2。注意2個高斯分布具有相同均方差Sigma，分別為Mu1,Mu2。

def ini_data(Sigma,Mu1,Mu2,k,N):

global X

global Mu

global Expectations

X = np.zeros((1,N))

Mu = np.random.random(2)

Expectations = np.zeros((N,k))

for i in xrange(0,N):

if np.random.random(1) > 0.5:

X[0,i] = np.random.normal()*Sigma + Mu1

else:

X[0,i] = np.random.normal()*Sigma + Mu2

if isdebug:

print "***********"

print u"初始觀測數據X："

print X

# EM算法：步驟1，計算E[zij]

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 Expectations

# EM算法：步驟2，求最大化E[zij]的參數Mu

def m_step(k,N):

global Expectations

global X

for j in xrange(0,k):

Numer = 0

Denom = 0

for i in xrange(0,N):

Numer += Expectations[i,j]*X[0,i]

Denom +=Expectations[i,j]

Mu[j] = Numer / Denom

# 算法迭代iter_num次，或達到精度Epsilon停止迭代

def run(Sigma,Mu1,Mu2,k,N,iter_num,Epsilon):

ini_data(Sigma,Mu1,Mu2,k,N)

print u"初始:", Mu

for i in range(iter_num):

Old_Mu = copy.deepcopy(Mu)

e_step(Sigma,k,N)

m_step(k,N)

print i,Mu

if sum(abs(Mu-Old_Mu)) < Epsilon:

break

if __name__ == '__main__':

run(6,40,20,2,1000,1000,0.0001)

plt.hist(X[0,:],50)

plt.show()

本代碼用于模擬k=2個正態分布的均值估計。其中ini_data(Sigma,Mu1,Mu2,k,N)函數用于生成訓練樣本，此訓練樣本時從兩個高斯分布中隨機生成的，其中高斯分布a均值Mu1=40、均方差Sigma=6，高斯分布b均值Mu2=20、均方差Sigma=6，生成的樣本分布如下圖所示。由于本問題中實現無法直接沖樣本數據中獲知兩個高斯分布參數，因此需要使用EM算法估算出具體Mu1、Mu2取值。

圖 1 ?樣本數據分布

在圖1的樣本數據下，在第11步時，迭代終止，EM估計結果為：

Mu=[ 40.55261688 ?19.34252468]

附：

總結

以上是生活随笔為你收集整理的python概率密度函数参数估计_EM算法求高斯混合模型参数估计——Python实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。