今天就跟大家聊聊有關(guān)使用python如何實現(xiàn)一個em算法,可能很多人都不太了解,為了讓大家更加了解,小編給大家總結(jié)了以下內(nèi)容,希望大家根據(jù)這篇文章可以有所收獲。
創(chuàng)新互聯(lián)服務(wù)項目包括斗門網(wǎng)站建設(shè)、斗門網(wǎng)站制作、斗門網(wǎng)頁制作以及斗門網(wǎng)絡(luò)營銷策劃等。多年來,我們專注于互聯(lián)網(wǎng)行業(yè),利用自身積累的技術(shù)優(yōu)勢、行業(yè)經(jīng)驗、深度合作伙伴關(guān)系等,向廣大中小型企業(yè)、政府機(jī)構(gòu)等提供互聯(lián)網(wǎng)行業(yè)的解決方案,斗門網(wǎng)站推廣取得了明顯的社會效益與經(jīng)濟(jì)效益。目前,我們服務(wù)的客戶以成都為中心已經(jīng)輻射到斗門省份的部分城市,未來相信會繼續(xù)擴(kuò)大服務(wù)區(qū)域并繼續(xù)獲得客戶的支持與信任!''' 數(shù)據(jù)集:偽造數(shù)據(jù)集(兩個高斯分布混合) 數(shù)據(jù)集長度:1000 ------------------------------ 運(yùn)行結(jié)果: ---------------------------- 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): ''' 初始化數(shù)據(jù)集 這里通過服從高斯分布的隨機(jī)函數(shù)來偽造數(shù)據(jù)集 :param mu0: 高斯0的均值 :param sigma0: 高斯0的方差 :param mu1: 高斯1的均值 :param sigma1: 高斯1的方差 :param alpha0: 高斯0的系數(shù) :param alpha1: 高斯1的系數(shù) :return: 混合了兩個高斯分布的數(shù)據(jù) ''' # 定義數(shù)據(jù)集長度為1000 length = 1000 # 初始化第一個高斯分布,生成數(shù)據(jù),數(shù)據(jù)長度為length * alpha系數(shù),以此來 # 滿足alpha的作用 data0 = np.random.normal(mu0, sigma0, int(length * alpha0)) # 第二個高斯分布的數(shù)據(jù) data1 = np.random.normal(mu1, sigma1, int(length * alpha1)) # 初始化總數(shù)據(jù)集 # 兩個高斯分布的數(shù)據(jù)混合后會放在該數(shù)據(jù)集中返回 dataSet = [] # 將第一個數(shù)據(jù)集的內(nèi)容添加進(jìn)去 dataSet.extend(data0) # 添加第二個數(shù)據(jù)集的數(shù)據(jù) dataSet.extend(data1) # 對總的數(shù)據(jù)集進(jìn)行打亂(其實不打亂也沒事,只不過打亂一下直觀上讓人感覺已經(jīng)混合了 # 讀者可以將下面這句話屏蔽以后看看效果是否有差別) random.shuffle(dataSet) #返回偽造好的數(shù)據(jù)集 return dataSet def calcGauss(dataSetArr, mu, sigmod): ''' 根據(jù)高斯密度函數(shù)計算值 依據(jù):“9.3.1 高斯混合模型” 式9.25 注:在公式中y是一個實數(shù),但是在EM算法中(見算法9.2的E步),需要對每個j 都求一次yjk,在本實例中有1000個可觀測數(shù)據(jù),因此需要計算1000次??紤]到 在E步時進(jìn)行1000次高斯計算,程序上比較不簡潔,因此這里的y是向量,在numpy 的exp中如果exp內(nèi)部值為向量,則對向量中每個值進(jìn)行exp,輸出仍是向量的形式。 所以使用向量的形式1次計算即可將所有計算結(jié)果得出,程序上較為簡潔 :param dataSetArr: 可觀測數(shù)據(jù)集 :param mu: 均值 :param sigmod: 方差 :return: 整個可觀測數(shù)據(jù)集的高斯分布密度(向量形式) ''' # 計算過程就是依據(jù)式9.25寫的,沒有別的花樣 result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2)) # 返回結(jié)果 return result def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1): ''' EM算法中的E步 依據(jù)當(dāng)前模型參數(shù),計算分模型k對觀數(shù)據(jù)y的響應(yīng)度 :param dataSetArr: 可觀測數(shù)據(jù)y :param alpha0: 高斯模型0的系數(shù) :param mu0: 高斯模型0的均值 :param sigmod0: 高斯模型0的方差 :param alpha1: 高斯模型1的系數(shù) :param mu1: 高斯模型1的均值 :param sigmod1: 高斯模型1的方差 :return: 兩個模型各自的響應(yīng)度 ''' # 計算y0的響應(yīng)度 # 先計算模型0的響應(yīng)度的分子 gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0) # 模型1響應(yīng)度的分子 gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1) # 兩者相加為E步中的分布 sum = gamma0 + gamma1 # 各自相除,得到兩個模型的響應(yīng)度 gamma0 = gamma0 / sum gamma1 = gamma1 / sum # 返回兩個模型響應(yīng)度 return gamma0, gamma1 def M_step(muo, mu1, gamma0, gamma1, dataSetArr): # 依據(jù)算法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): ''' 根據(jù)EM算法進(jìn)行參數(shù)估計 算法依據(jù)“9.3.2 高斯混合模型參數(shù)估計的EM算法” 算法9.2 :param dataSetList:數(shù)據(jù)集(可觀測數(shù)據(jù)) :param iter: 迭代次數(shù) :return: 估計的參數(shù) ''' # 將可觀測數(shù)據(jù)y轉(zhuǎn)換為數(shù)組形式,主要是為了方便后續(xù)運(yùn)算 dataSetArr = np.array(dataSetList) # 步驟1:對參數(shù)取初值,開始迭代 alpha0 = 0.5 mu0 = 0 sigmod0 = 1 alpha1 = 0.5 mu1 = 1 sigmod1 = 1 # 開始迭代 step = 0 while (step < iter): # 每次進(jìn)入一次迭代后迭代次數(shù)加1 step += 1 # 步驟2:E步:依據(jù)當(dāng)前模型參數(shù),計算分模型k對觀測數(shù)據(jù)y的響應(yīng)度 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) # 迭代結(jié)束后將更新后的各參數(shù)返回 return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 if __name__ == '__main__': start = time.time() # 設(shè)置兩個高斯模型進(jìn)行混合,這里是初始化兩個模型各自的參數(shù) # 見“9.3 EM算法在高斯混合模型學(xué)習(xí)中的應(yīng)用” # alpha是“9.3.1 高斯混合模型” 定義9.2中的系數(shù)α # mu0是均值μ # sigmod是方差σ # 在設(shè)置上兩個alpha的和必須為1,其他沒有什么具體要求,符合高斯定義就可以 alpha0 = 0.3 # 系數(shù)α mu0 = -2 # 均值μ sigmod0 = 0.5 # 方差σ alpha1 = 0.7 # 系數(shù)α mu1 = 0.5 # 均值μ sigmod1 = 1 # 方差σ # 初始化數(shù)據(jù)集 dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1) #打印設(shè)置的參數(shù) 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算法,進(jìn)行參數(shù)估計 alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList) # 打印參數(shù)預(yù)測結(jié)果 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)