【machine learning】GMM演算法（Python版）
本文參考CSDN大神的博文，並在講述中引入自己的理解，純粹理清思路，並將程式碼改為了Python版本。（在更改的過程中，一方面理清自己對GMM的理解，一方面學習了numpy的應用，不過也許是Python粉指數超標才覺得有必要改(⊙o⊙)）
一、GMM模型
事實上，GMM 和 k-means 很像，不過 GMM 是學習出一些概率密度函式來（所以 GMM 除了用在 clustering 上之外，還經常被用於 density estimation ），簡單地說，k-means 的結果是每個資料點被 assign 到其中某一個 cluster 了，而 GMM 則給出這些資料點被 assign 到每個 cluster 的概率，又稱作 soft assignment 。
得出一個概率有很多好處，因為它的資訊量比簡單的一個結果要多，比如，我可以把這個概率轉換為一個 score ，表示演算法對自己得出的這個結果的把握，參考pluskid大神博文
- 也許我可以對同一個任務，用多個方法得到結果，最後選取“把握”最大的那個結果；
- 另一個很常見的方法是在諸如疾病診斷之類的場所，機器對於那些很容易分辨的情況（患病或者不患病的概率很高）可以自動區分，而對於那種很難分辨的情況，比如，49% 的概率患病，51% 的概率正常，如果僅僅簡單地使用 50% 的閾值將患者診斷為“正常”的話，風險是非常大的，因此，在機器對自己的結果把握很小的情況下，會“拒絕發表評論”，而把這個任務留給有經驗的醫生去解決。
每個GMM由K個Gaussian分佈組成，每個Gaussian稱為一個“Component”，這些Component 線性加成在一起就組成了GMM 的概率密度函式：
根據上面的式子，如果我們要從 GMM 的分佈中隨機地取一個點的話，實際上可以分為兩步：首先隨機地在這 K個Gaussian Component 之中選一個，每個 Component 被選中的概率實際上就是它的係數 pi(k) ，選中了 Component 之後，再單獨地考慮從這個 Component 的分佈中選取一個點就可以了──這裡已經回到了普通的 Gaussian 分佈，轉化為了已知的問題。
那麼如何用 GMM 來做 clustering 呢？其實很簡單，現在我們有了資料，假定它們是由 GMM 生成出來的，那麼我們只要根據資料推出 GMM 的概率分佈來就可以了，然後 GMM 的 K 個 Component 實際上就對應了 K 個 cluster 了。根據資料來推算概率密度通常被稱作 density estimation ，特別地，當我們在已知（或假定）了概率密度函式的形式，而要估計其中的引數的過程被稱作“引數估計”。
二、引數與似然函式
現在假設我們有 N 個數據點，並假設它們服從某個分佈（記作 p(x)），現在要確定裡面的一些引數的值，例如，在 GMM 中，我們就需要確定 影響因子pi(k)、各類均值pMiu(k) 和 各類協方差pSigma(k) 這些引數。 我們的想法是，找到這樣一組引數，它所確定的概率分佈生成這些給定的資料點的概率最大，而這個概率實際上就等於 ，我們把這個乘積稱作似然函式 (Likelihood Function)。通常單個點的概率都很小，許多很小的數字相乘起來在計算機裡很容易造成浮點數下溢，因此我們通常會對其取對數，把乘積變成加和
下面讓我們來看一看 GMM 的 log-likelihood function ：
由於在對數函式裡面又有加和，我們沒法直接用求導解方程的辦法直接求得最大值。為了解決這個問題，我們採取之前從 GMM 中隨機選點的辦法：分成兩步，實際上也就類似於K-means 的兩步。
三、演算法流程
1、估計資料由每個 Component 生成的概率（並不是每個 Component 被選中的概率）：對於每個資料
其中N（xi|μk,Σk）就是後驗概率
2、通過極大似然估計可以通過求到令引數=0得到引數pMiu，pSigma的值。
其中
3、重複迭代前面兩步，直到似然函式的值收斂為止。
宣告：這裡完全可以對照EM演算法的兩步走，對應關係如下：
- 均值和方差對應θ
- Component 生成的概率為隱藏變數
- 最大似然函式L(θ)=
四、GMM聚類程式碼
GMM.py
#! /usr/bin/env python
#coding=utf-8
from numpy import *
import pylab
import random,math
def loadDataSet(fileName): #general function to parse tab -delimited floats
dataMat = [] #assume last column is target value
fr = open(fileName)
for line in fr.readlines():
curLine = line.strip().split('\t')
fltLine = map(float,curLine) #map all elements to float()
dataMat.append(fltLine)
return dataMat
def gmm(file, K_or_centroids):
# ============================================================
# Expectation-Maximization iteration implementation of
# Gaussian Mixture Model.
#
# PX = GMM(X, K_OR_CENTROIDS)
# [PX MODEL] = GMM(X, K_OR_CENTROIDS)
#
# - X: N-by-D data matrix.
# - K_OR_CENTROIDS: either K indicating the number of
# components or a K-by-D matrix indicating the
# choosing of the initial K centroids.
#
# - PX: N-by-K matrix indicating the probability of each
# component generating each point.
# - MODEL: a structure containing the parameters for a GMM:
# MODEL.Miu: a K-by-D matrix.
# MODEL.Sigma: a D-by-D-by-K matrix.
# MODEL.Pi: a 1-by-K vector.
# ============================================================
## Generate Initial Centroids
threshold = 1e-15
dataMat = mat(loadDataSet(file))
[N, D] = shape(dataMat)
K_or_centroids = 2
# K_or_centroids可以是一個整數，也可以是k個質心的二維列向量
if shape(K_or_centroids)==(): #if K_or_centroid is a 1*1 number
K = K_or_centroids
Rn_index = range(N)
random.shuffle(Rn_index) #random index N samples
centroids = dataMat[Rn_index[0:K], :]; #generate K random centroid
else: # K_or_centroid is a initial K centroid
K = size(K_or_centroids)[0];
centroids = K_or_centroids;
## initial values
[pMiu,pPi,pSigma] = init_params(dataMat,centroids,K,N,D)
Lprev = -inf #上一次聚類的誤差
# EM Algorithm
while True:
# Estimation Step
Px = calc_prob(pMiu,pSigma,dataMat,K,N,D)
# new value for pGamma(N*k), pGamma(i,k) = Xi由第k個Gaussian生成的概率
# 或者說xi中有pGamma(i,k)是由第k個Gaussian生成的
pGamma = mat(array(Px) * array(tile(pPi, (N, 1)))) #分子 = pi(k) * N(xi | pMiu(k), pSigma(k))
pGamma = pGamma / tile(sum(pGamma, 1), (1, K)) #分母 = pi(j) * N(xi | pMiu(j), pSigma(j))對所有j求和
## Maximization Step - through Maximize likelihood Estimation
#print 'dtypeddddddddd:',pGamma.dtype
Nk = sum(pGamma, 0) #Nk(1*k) = 第k個高斯生成每個樣本的概率的和，所有Nk的總和為N。
# update pMiu
pMiu = mat(diag((1/Nk).tolist()[0])) * (pGamma.T) * dataMat #update pMiu through MLE(通過令導數 = 0得到)
pPi = Nk/N
# update k個 pSigma
print 'kk=',K
for kk in range(K):
Xshift = dataMat-tile(pMiu[kk], (N, 1))
Xshift.T * mat(diag(pGamma[:, kk].T.tolist()[0])) * Xshift / 2
pSigma[:, :, kk] = (Xshift.T * \
mat(diag(pGamma[:, kk].T.tolist()[0])) * Xshift) / Nk[kk]
# check for convergence
L = sum(log(Px*(pPi.T)))
if L-Lprev < threshold:
break
Lprev = L
return Px
def init_params(X,centroids,K,N,D):
pMiu = centroids #k*D, 即k類的中心點
pPi = zeros([1, K]) #k類GMM所佔權重（influence factor）
pSigma = zeros([D, D, K]) #k類GMM的協方差矩陣，每個是D*D的
# 距離矩陣，計算N*K的矩陣（x-pMiu）^2 = x^2+pMiu^2-2*x*Miu
#x^2, N*1的矩陣replicateK列\#pMiu^2，1*K的矩陣replicateN行
distmat = tile(sum(power(X,2), 1),(1, K)) + \
tile(transpose(sum(power(pMiu,2), 1)),(N, 1)) - \
2*X*transpose(pMiu)
labels = distmat.argmin(1) #Return the minimum from each row
# 獲取k類的pPi和協方差矩陣
for k in range(K):
boolList = (labels==k).tolist()
indexList = [boolList.index(i) for i in boolList if i==[True]]
Xk = X[indexList, :]
#print cov(Xk)
# 也可以用shape(XK)[0]
pPi[0][k] = float(size(Xk, 0))/N
pSigma[:, :, k] = cov(transpose(Xk))
return pMiu,pPi,pSigma
# 計算每個資料由第k類生成的概率矩陣Px
def calc_prob(pMiu,pSigma,X,K,N,D):
# Gaussian posterior probability
# N(x|pMiu,pSigma) = 1/((2pi)^(D/2))*(1/(abs(sigma))^0.5)*exp(-1/2*(x-pMiu)'pSigma^(-1)*(x-pMiu))
Px = mat(zeros([N, K]))
for k in range(K):
Xshift = X-tile(pMiu[k, :],(N, 1)) #X-pMiu
#inv_pSigma = mat(pSigma[:, :, k]).I
inv_pSigma = linalg.pinv(mat(pSigma[:, :, k]))
tmp = sum(array((Xshift*inv_pSigma)) * array(Xshift), 1) # 這裡應變為一列數
tmp = mat(tmp).T
#print linalg.det(inv_pSigma),'54545'
Sigema = linalg.det(mat(inv_pSigma))
if Sigema < 0:
Sigema=0
coef = power((2*(math.pi)),(-D/2)) * sqrt(Sigema)
Px[:, k] = coef * exp(-0.5*tmp)
return Pxmain.py
#! /usr/bin/env python
#coding=utf-8
import GMM
'''
def showFigure(dataMat,k,clusterAssment):
tag=['go','or','yo','ko']
for i in range(k):
datalist = dataMat[nonzero(clusterAssment[:,0].A==i)[0]]
pylab.plot(datalist[:,0],datalist[:,1],tag[i])
pylab.show()
'''
if __name__ == '__main__':
GMM.gmm('testSet.txt',2)說明：