'''
數據集:偽造數據集(兩個高斯分布混合)
數據集長度:1000
------------------------------
運行結果:
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''
import numpy as np
import random
import math
import time
def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
'''
初始化數據集
這里通過服從高斯分布的隨機函數來偽造數據集
:param mu0: 高斯0的均值
:param sigma0: 高斯0的方差
:param mu1: 高斯1的均值
:param sigma1: 高斯1的方差
:param alpha0: 高斯0的系數
:param alpha1: 高斯1的系數
:return: 混合了兩個高斯分布的數據
'''
# 定義數據集長度為1000
length = 1000
# 初始化第一個高斯分布,生成數據,數據長度為length * alpha系數,以此來
# 滿足alpha的作用
data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
# 第二個高斯分布的數據
data1 = np.random.normal(mu1, sigma1, int(length * alpha1))
# 初始化總數據集
# 兩個高斯分布的數據混合后會放在該數據集中返回
dataSet = []
# 將第一個數據集的內容添加進去
dataSet.extend(data0)
# 添加第二個數據集的數據
dataSet.extend(data1)
# 對總的數據集進行打亂(其實不打亂也沒事,只不過打亂一下直觀上讓人感覺已經混合了
# 讀者可以將下面這句話屏蔽以后看看效果是否有差別)
random.shuffle(dataSet)
#返回偽造好的數據集
return dataSet
def calcGauss(dataSetArr, mu, sigmod):
'''
根據高斯密度函數計算值
依據:“9.3.1 高斯混合模型” 式9.25
注:在公式中y是一個實數,但是在EM算法中(見算法9.2的E步),需要對每個j
都求一次yjk,在本實例中有1000個可觀測數據,因此需要計算1000次。考慮到
在E步時進行1000次高斯計算,程序上比較不簡潔,因此這里的y是向量,在numpy
的exp中如果exp內部值為向量,則對向量中每個值進行exp,輸出仍是向量的形式。
所以使用向量的形式1次計算即可將所有計算結果得出,程序上較為簡潔
:param dataSetArr: 可觀測數據集
:param mu: 均值
:param sigmod: 方差
:return: 整個可觀測數據集的高斯分布密度(向量形式)
'''
# 計算過程就是依據式9.25寫的,沒有別的花樣
result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
# 返回結果
return result
def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
'''
EM算法中的E步
依據當前模型參數,計算分模型k對觀數據y的響應度
:param dataSetArr: 可觀測數據y
:param alpha0: 高斯模型0的系數
:param mu0: 高斯模型0的均值
:param sigmod0: 高斯模型0的方差
:param alpha1: 高斯模型1的系數
:param mu1: 高斯模型1的均值
:param sigmod1: 高斯模型1的方差
:return: 兩個模型各自的響應度
'''
# 計算y0的響應度
# 先計算模型0的響應度的分子
gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
# 模型1響應度的分子
gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)
# 兩者相加為E步中的分布
sum = gamma0 + gamma1
# 各自相除,得到兩個模型的響應度
gamma0 = gamma0 / sum
gamma1 = gamma1 / sum
# 返回兩個模型響應度
return gamma0, gamma1
def M_step(muo, mu1, gamma0, gamma1, dataSetArr):
# 依據算法9.2計算各個值
# 這里沒什么花樣,對照書本公式看看這里就好了
mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)
sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))
sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))
alpha0_new = np.sum(gamma0) / len(gamma0)
alpha1_new = np.sum(gamma1) / len(gamma1)
# 將更新的值返回
return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new
def EM_Train(dataSetList, iter=500):
'''
根據EM算法進行參數估計
算法依據“9.3.2 高斯混合模型參數估計的EM算法” 算法9.2
:param dataSetList:數據集(可觀測數據)
:param iter: 迭代次數
:return: 估計的參數
'''
# 將可觀測數據y轉換為數組形式,主要是為了方便后續運算
dataSetArr = np.array(dataSetList)
# 步驟1:對參數取初值,開始迭代
alpha0 = 0.5
mu0 = 0
sigmod0 = 1
alpha1 = 0.5
mu1 = 1
sigmod1 = 1
# 開始迭代
step = 0
while (step < iter):
# 每次進入一次迭代后迭代次數加1
step += 1
# 步驟2:E步:依據當前模型參數,計算分模型k對觀測數據y的響應度
gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
# 步驟3:M步
mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
# 迭代結束后將更新后的各參數返回
return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1
if __name__ == '__main__':
start = time.time()
# 設置兩個高斯模型進行混合,這里是初始化兩個模型各自的參數
# 見“9.3 EM算法在高斯混合模型學習中的應用”
# alpha是“9.3.1 高斯混合模型” 定義9.2中的系數α
# mu0是均值μ
# sigmod是方差σ
# 在設置上兩個alpha的和必須為1,其他沒有什么具體要求,符合高斯定義就可以
alpha0 = 0.3 # 系數α
mu0 = -2 # 均值μ
sigmod0 = 0.5 # 方差σ
alpha1 = 0.7 # 系數α
mu1 = 0.5 # 均值μ
sigmod1 = 1 # 方差σ
# 初始化數據集
dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)
#打印設置的參數
print('---------------------------')
print('the Parameters set is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
# 開始EM算法,進行參數估計
alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)
# 打印參數預測結果
print('----------------------------')
print('the Parameters predict is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
# 打印時間
print('----------------------------')
print('time span:', time.time() - start)