【machine learning】GMM算法（Python版）
一、GMM模型
事實上，GMM 和 k-means 很像，不過 GMM 是學(xué)習(xí)出一些概率密度函數(shù)來（所以 GMM 除了用在 clustering 上之外，還經(jīng)常被用于 density estimation ），簡單地說，k-means 的結(jié)果是每個數(shù)據(jù)點(diǎn)被 assign 到其中某一個 cluster 了，而 GMM 則給出這些數(shù)據(jù)點(diǎn)被 assign 到每個 cluster 的概率，又稱作 soft assignment 。
得出一個概率有很多好處，因為它的信息量比簡單的一個結(jié)果要多，比如，我可以把這個概率轉(zhuǎn)換為一個 score ，表示算法對自己得出的這個結(jié)果的把握，參考pluskid大神博文
- 也許我可以對同一個任務(wù)，用多個方法得到結(jié)果，最后選取“把握”最大的那個結(jié)果；
- 另一個很常見的方法是在諸如疾病診斷之類的場所，機(jī)器對于那些很容易分辨的情況（患病或者不患病的概率很高）可以自動區(qū)分，而對于那種很難分辨的情況，比如，49% 的概率患病，51% 的概率正常，如果僅僅簡單地使用 50% 的閾值將患者診斷為“正常”的話，風(fēng)險是非常大的，因此，在機(jī)器對自己的結(jié)果把握很小的情況下，會“拒絕發(fā)表評論”，而把這個任務(wù)留給有經(jīng)驗的醫(yī)生去解決。
每個GMM由K個Gaussian分布組成，每個Gaussian稱為一個“Component”，這些Component 線性加成在一起就組成了GMM 的概率密度函數(shù)：
根據(jù)上面的式子，如果我們要從 GMM 的分布中隨機(jī)地取一個點(diǎn)的話，實際上可以分為兩步：首先隨機(jī)地在這 K個Gaussian Component 之中選一個，每個 Component 被選中的概率實際上就是它的系數(shù) pi(k) ，選中了 Component 之后，再單獨(dú)地考慮從這個 Component 的分布中選取一個點(diǎn)就可以了──這里已經(jīng)回到了普通的 Gaussian 分布，轉(zhuǎn)化為了已知的問題。
那么如何用 GMM 來做 clustering 呢？其實很簡單，現(xiàn)在我們有了數(shù)據(jù)，假定它們是由 GMM 生成出來的，那么我們只要根據(jù)數(shù)據(jù)推出 GMM 的概率分布來就可以了，然后 GMM 的 K 個 Component 實際上就對應(yīng)了 K 個 cluster 了。根據(jù)數(shù)據(jù)來推算概率密度通常被稱作 density estimation ，特別地，當(dāng)我們在已知（或假定）了概率密度函數(shù)的形式，而要估計其中的參數(shù)的過程被稱作“參數(shù)估計”。
二、參數(shù)與似然函數(shù)
現(xiàn)在假設(shè)我們有 N 個數(shù)據(jù)點(diǎn)，并假設(shè)它們服從某個分布（記作 p(x)），現(xiàn)在要確定里面的一些參數(shù)的值，例如，在 GMM 中，我們就需要確定 影響因子pi(k)、各類均值pMiu(k) 和 各類協(xié)方差pSigma(k) 這些參數(shù)。 我們的想法是，找到這樣一組參數(shù)，它所確定的概率分布生成這些給定的數(shù)據(jù)點(diǎn)的概率最大，而這個概率實際上就等于 ，我們把這個乘積稱作似然函數(shù) (Likelihood Function)。通常單個點(diǎn)的概率都很小，許多很小的數(shù)字相乘起來在計算機(jī)里很容易造成浮點(diǎn)數(shù)下溢，因此我們通常會對其取對數(shù)，把乘積變成加和 ∑ N i = 1 log p ( x i ) ∑i=1Nlog?p(xi) ，得到log-likelihood function 。接下來我們只要將這個函數(shù)最大化（通常的做法是求導(dǎo)并令導(dǎo)數(shù)等于零，然后解方程），亦即找到這樣一組參數(shù)值，它讓似然函數(shù)取得最大值，我們就認(rèn)為這是最合適的參數(shù)，這樣就完成了參數(shù)估計的過程。
下面讓我們來看一看 GMM 的 log-likelihood function ：
由于在對數(shù)函數(shù)里面又有加和，我們沒法直接用求導(dǎo)解方程的辦法直接求得最大值。為了解決這個問題，我們采取之前從 GMM 中隨機(jī)選點(diǎn)的辦法：分成兩步，實際上也就類似于K-means 的兩步。
三、算法流程
1、估計數(shù)據(jù)由每個 Component 生成的概率（并不是每個 Component 被選中的概率）：對于每個數(shù)據(jù) x i xi 來說，它由第 k 個 Component 生成的概率為
其中N（xi|μk,Σk）就是后驗概率
2、通過極大似然估計可以通過求到令參數(shù)=0得到參數(shù)pMiu，pSigma的值。
其中 N k = ∑ N i = 1 γ ( i , k ) Nk=∑i=1Nγ(i,k) (這里的順理成章要考證起來有很多要說)
3、重復(fù)迭代前面兩步，直到似然函數(shù)的值收斂為止。
聲明：這里完全可以對照EM算法的兩步走，對應(yīng)關(guān)系如下：
- 均值和方差對應(yīng)θ
- Component 生成的概率為隱藏變量
-
最大似然函數(shù)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可以是一個整數(shù)，也可以是k個質(zhì)心的二維列向量
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(通過令導(dǎo)數(shù) = 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類的中心點(diǎn)
pPi = zeros([
1
, K])
#k類GMM所占權(quán)重（influence factor）
pSigma = zeros([D, D, K])
#k類GMM的協(xié)方差矩陣，每個是D*D的
# 距離矩陣，計算N*K的矩陣（x-pMiu）^2 = x^2+pMiu^2-2*x*Miu
#x^2, N*1的矩陣replicateK列\(zhòng)#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和協(xié)方差矩陣
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
# 計算每個數(shù)據(jù)由第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
)
# 這里應(yīng)變?yōu)橐涣袛?shù)
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
Px
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
main.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
)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
說明：
-
程序由matlab代碼改過了，利用了numpy庫的函數(shù)，需要下載軟件可參考鏈接
-
其中的分析數(shù)據(jù)源我用了KMeans算法中的數(shù)據(jù)，不過遇到了奇異矩陣的問題（因為數(shù)據(jù)分布本來就有悖于高斯分布），即使用了偽逆也沒解決問題，看之后學(xué)習(xí)能否回頭再解決。也可查看博文評論圍觀別人的解決方法
-
matlab在矩陣的處理上的確優(yōu)于Python太多了，一方面不用導(dǎo)入庫，也沒有存在array,mat等類的轉(zhuǎn)換和眾多函數(shù)的比較和考慮。不過Python綜合應(yīng)用多，上至Web下至PC，從測試到開發(fā)都有很多人用，相比matlab還是做測試多一點(diǎn)。
更多文章、技術(shù)交流、商務(wù)合作、聯(lián)系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯(lián)系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點(diǎn)擊下面給點(diǎn)支持吧，站長非常感激您！手機(jī)微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點(diǎn)擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
