Python：EM(期望極大算法)實(shí)戰(zhàn)

　　先來詳細(xì)的講一下EM算法。

　　前提準(zhǔn)備

　　Jupyter notebook 或 Pycharm

　　火狐瀏覽器或Chrome瀏覽器

　　win7或win10電腦一臺

　　網(wǎng)盤提取csv數(shù)據(jù)

　　需求分析

　　實(shí)現(xiàn)高斯混合模型的 EM 算法(GMM_EM)

　　高斯混合模型是多個高斯模型的線性疊加而成的，高斯混合模型的概率分布表示如下：

　　其中，k表示模型的個數(shù)，αkα_kαk 是第 k 個模型的系數(shù)，表示出現(xiàn)該模型的概率，?(x;μk,Σk) 是第 k 個高斯模型的概率分布。

　　E步：樣本 xix_ixi來自于第 k 個模型的概率，我們把這個概率稱為模型 k 對樣本 xix_ixi 的“責(zé)任”，也叫“響應(yīng)度”，記作 γ(ik)γ_(ik)γ(ik)，計(jì)算公式如下：

　　M步：根據(jù)樣本和當(dāng)前 γ 矩陣重新估計(jì)參數(shù)，注意這里 x 為列向量

　　【目標(biāo)】給定一堆沒有標(biāo)簽的樣本和模型個數(shù) K，以此求得混合模型的參數(shù)，然后就可以用這個模型來對樣本進(jìn)行聚類。

　　python代碼如下：

　　import numpy as np

　　import matplotlib.pyplot as plt

　　from scipy.stats import multivariate_normal #本問題考慮的是高斯混合模型，所以導(dǎo)入多元高斯分布multivariate_normal

　　def prob_Y_k(Y,mu_k,cov_k): #Y為樣本矩陣

　　norm = multivariate_normal(mean = mu_k , cov = cov_k) #生成多元正太分布，mu為第k個模型的均值，cov為第k個模型的協(xié)方差矩陣(協(xié)方差矩陣必須是實(shí)對稱矩陣)

　　return norm.pdf(Y) #返回樣本Y來自于第k個模型的概率

　　def Estep(Y,mu,cov,alpha): #Y為樣本矩陣，alpha為權(quán)重

　　N = Y.shape[0] #樣本數(shù)

　　K = alpha.shape[0] #模型數(shù)

　　assert N>1 , "There must be more than one sample!" #為避免單個樣本導(dǎo)致返回的結(jié)果的類型不一致，因此要求樣本數(shù)必須大于一

　　assert K>1 , "There must be more than one gaussian model!" #為避免單個模型結(jié)果的類型不一致，因此要求模型須大于一

　　gamma = np.mat(np.zeros((N,K))) #初始化響應(yīng)度矩陣，行對應(yīng)樣本數(shù)，列對應(yīng)模型數(shù)

　　prob = np.zeros((N,K)) #初始化所有樣本出現(xiàn)的概率矩陣，行對應(yīng)樣本數(shù)，列對應(yīng)響應(yīng)度

　　for k in range(K):

　　prob[:,k] = prob_Y_k(Y,mu[k],cov[k]) #第k個模型的概率prob_Y_k

　　prob = np.mat(prob) #K個prob放入數(shù)組中

　　for k in range(K):

　　gamma[:,k] = alpha[k] * prob[:,k] #計(jì)算模型k對樣本i的響應(yīng)度

　　for i in range(N):

　　gamma[i,:] /= np.sum(gamma[i,:]) #第i個樣本的占總樣本的響應(yīng)程度

　　return gamma #gamma為響應(yīng)度矩陣

　　def Mstep(Y,gamma): #傳入樣本矩陣Y和Estep得到的gamma響應(yīng)度矩陣

　　N, D = Y.shape #N為樣本數(shù)，D為特征數(shù)

　　K = gamma.shape[1] #模型數(shù)

　　mu = np.zeros((K,D)) #初始化參數(shù)均值mu，每個模型的D維各有均值故mu的矩陣為K行D列

　　cov = [] #初始化參數(shù)協(xié)方差矩陣

　　alpha = np.zeros(K) # 初始化權(quán)重數(shù)組，每個模型都有權(quán)值

　　#接下來是更新每個模型的參數(shù)

　　for k in range(K):

　　Nk = np.sum(gamma[:,k]) #第k個模型所有樣本的響應(yīng)度之和

　　mu[k,:] = np.sum(np.multiply(Y, gamma[:,k]),axis=0)/Nk #更新參數(shù)均值mu，對每個特征求均值

　　cov_k = (Y - mu[k]).T * np.multiply((Y - mu[k]), gamma[:,k]) / Nk #更新cov

　　cov = np.append(cov_k)

　　alpha[k] = Nk / N

　　cov = np.array(cov)

　　return mu, cov, alpha

　　def normalize_data(Y): #將所有數(shù)據(jù)進(jìn)行歸一化處理，

　　for i in range(Y.shape[1]):

　　max_data = Y[:,i].max()

　　min_data = Y[:,i].min()

　　Y[:,i] = (Y[:,i] - min_data)/(max_data - min_data) #此處用到min-max歸一化

　　debug("Data Normalized")

　　return Y

　　def init_params(shape,K): #在執(zhí)行該算法之前，需要先給出一個初始化的模型參數(shù)。我們讓每個模型的μ為隨機(jī)值,Σ 為單位矩陣，α 為 1/K，即每個模型初始時都是等概率出現(xiàn)的。

　　N, D = shape

　　mu = np.random.rand(K, D) #生成一個K行D列的[0,1)之間的數(shù)組

　　cov = np.array([np.eye(D)] * K) #生成K個D維的對角矩陣

　　alpha = np.array([1.0 / K] * K) #生成K個權(quán)重

　　debug("Parameters initialized.")

　　debug("mu:",mu, "cov:",cov ,"alpha:",alpha,sep = "\n" )

　　return mu, cov, alpha

　　def GMM_EM(Y, K, times): #高斯混合EM算法,Y為給定樣本矩陣，K為模型個數(shù)，times為迭代次數(shù)，目的是求該模型的參數(shù)

　　Y = normalize_data(Y) #調(diào)用前面定義的normalize_data函數(shù)，歸一化樣本矩陣Y

　　mu, cov, alpha = init_params(Y.shape, K) #調(diào)用init_params函數(shù)得到初始化的參數(shù)mu,cov,alpha

　　for i in range(times):鄭州婦科醫(yī)院 http://m.zyfuke.com/

　　gamma = Estep(Y, mu, cov, alpha) #調(diào)用Estep得到響應(yīng)度矩陣

　　mu, cov, alpha = Mstep(Y, gamma) #調(diào)用Mstep得到更新后的參數(shù)mu,cov,alpha

　　debug("{sep} Result {sep}".format(sep="-"*20))

　　debug("mu:", mu , "cov:",cov , "alpha:",alpha , sep="\n")

　　return mu,cov,alpha

　　import matplotlib.pyplot as plt

　　from gmm import *

　　DEBUG = True

　　Y = np.loadtxt("gmm.data") #載入數(shù)據(jù)

　　matY = np.matrix(Y ,copy = True)

　　K = 2 #模型個數(shù)(相當(dāng)于聚類的類別個數(shù))

　　mu, cov, alpha = GMM_EM(matY , K , 100) #調(diào)用GMM_EM函數(shù)，計(jì)算GMM模型參數(shù)

　　N = Y.shape[0]

　　gamma = Estep(matY, mu, cov, alpha) #求當(dāng)前模型參數(shù)下，各模型對樣本的響應(yīng)矩陣

　　category = gamma.argmax(axis = 1).flatten().tolist()[0] #對每個樣本，求響應(yīng)度最大的模型下標(biāo)，作為其類別標(biāo)識

　　class1 = np.array([Y[i] for i in range(N) if category[i] == 0]) #將每個樣本放入對應(yīng)樣本的列表中

　　class2 = np.array([Y[i] for i in range(N) if category[i] == 1])

　　plt.plot(class1[:,0],class1[:,1], 'rs' ,label = "class1")

　　plt.plot(class2[:,0],class2[:,1], 'bo' ,label = "class2")

　　plt.legend(loc = "best")

　　plt.title("GMM Clustering By EM Algorithm")

　　plt.show()

　　import numpy as np

　　import matplotlib.pyplot as plt

　　cov1 = np.mat("0.3 0 ; 0 0.1") #2維協(xié)方差矩陣(必須是對角矩陣)

　　cov2 = np.mat("0.2 0 ; 0 0.3")

　　mu1 = np.array([0,1])

　　mu2 = np.array([2,1])

　　sample = np.zeros((100,2)) #初始化100個樣本，樣本特征為2

　　sample[:30, :] = np.random.multivariate_normal(mean=mu1, cov=cov1, size=30) #生成多元正態(tài)分布矩陣

　　sample[30:, :] = np.random.multivariate_normal(mean=mu2, cov=cov2, size=70)

　　np.savetxt("sample.data",sample) # 將array保存到txt文件中

　　plt.plot(sample[:30, 0], sample[:30, 1], "bo") #30個樣本用藍(lán)色圓圈標(biāo)記

　　plt.plot(sample[30:, 0], sample[30:, 1], "rs") #70個樣本用紅色方塊標(biāo)記

　　plt.title("sample_data")

　　plt.show()