作者：金良（golden1314521@gmail.com） csdn博客： http://blog.csdn.net/u012176591

需要整理后的代碼文件和數(shù)據(jù)請(qǐng)移步 http://download.csdn.net/detail/u012176591/8748673

1.高斯分布公式及圖像示例

定義在D-維連續(xù)空間的高斯分布概率密度的表達(dá)式
N(x|μ,Σ)=1(2π)D/21|Σ|1/2exp{?12(x?μ)TΣ?1(x?μ)}

其等高線所形成的形狀與協(xié)方差矩陣Σ 密切相關(guān)，如下所示，后面的代碼中有各個(gè)圖像的對(duì)應(yīng)的高斯分布的參數(shù)。

2. 高斯分布概率密度熱力圖

代碼如下：

fig,axes = plt.subplots(nrows=3,ncols=1,figsize=(4,12))

# 標(biāo)準(zhǔn)圓形
mean = [0,0]
cov = [[1,0],
       [0,1]] 
x,y = np.random.multivariate_normal(mean,cov,5000).T
axes[0].plot(x,y,'x')
axes[0].set_xlim(-6,6) 
axes[0].set_ylim(-6,6) 

# 橢圓，橢圓的軸向與坐標(biāo)平行
mean = [0,0]
cov = [[0.5,0],
       [0,3]] 
x,y = np.random.multivariate_normal(mean,cov,5000).T
axes[1].plot(x,y,'x')
axes[1].set_xlim(-6,6) 
axes[1].set_ylim(-6,6) 

# 橢圓，但是橢圓的軸與坐標(biāo)軸不一定平行
mean = [0,0]
cov = [[1,2.3],
       [2.3,1.4]] 
x,y = np.random.multivariate_normal(mean,cov,5000).T
axes[2].plot(x,y,'x'); plt.axis('equal')
axes[2].set_xlim(-6,6) 
axes[2].set_ylim(-6,6) 
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
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

我們?cè)谙旅娴母咚够旌夏Ｐ椭胁捎糜玫谌N協(xié)方差矩陣，即概率密度的等高線是橢圓，且軸向不一定與坐標(biāo)軸平行。

下圖是高斯密度函數(shù)的熱圖：

以下是作圖代碼

# 自定義的高維高斯分布概率密度函數(shù)
def gaussian(x,mean,cov):    
    dim = np.shape(cov)[0] #維度
    covdet = np.linalg.det(cov+np.eye(dim)*0.01) #協(xié)方差矩陣的秩
    covinv = np.linalg.inv(cov+np.eye(dim)*0.01) #協(xié)方差矩陣的逆
    xdiff = x - mean
    #概率密度
    prob = 1.0/np.power(2*np.pi,1.0*2/2)/np.sqrt(np.abs(covdet))*np.exp(-1.0/2*np.dot(np.dot(xdiff,covinv),xdiff))

    return prob


#作二維高斯概率密度函數(shù)的熱力圖
mean = [0,0]
cov = [[1,2.3],
       [2.3,1.4]] 
x,y = np.random.multivariate_normal(mean,cov,5000).T
cov = np.cov(x,y) #由真實(shí)數(shù)據(jù)計(jì)算得到的協(xié)方差矩陣，而不是自己任意設(shè)定
n=200
x = np.linspace(-6,6,n)
y = np.linspace(-6,6,n)
xx,yy = np.meshgrid(x, y)
zz = np.zeros((n,n))
for i in range(n):
    for j in range(n):
        zz[i][j] = gaussian(np.array([xx[i][j],yy[i][j]]),mean,cov)

gci = plt.imshow(zz,origin='lower') # 選項(xiàng)origin='lower' 防止tuixan圖像顛倒
plt.xticks([5,100,195],[-5,0,5])
plt.yticks([5,100,195],[-5,0,5])
plt.title(u'高斯函數(shù)的熱力圖',{'fontname':'STFangsong','fontsize':18})
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
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

３.高斯混合模型實(shí)現(xiàn)代碼：

下面是幾個(gè)功能函數(shù)，在主函數(shù)中被調(diào)用

# 計(jì)算概率密度，
# 參數(shù)皆為array類型，過(guò)程中參數(shù)不變
def gaussian(x,mean,cov):

    dim = np.shape(cov)[0] #維度
    #之所以加入單位矩陣是為了防止行列式為0的情況
    covdet = np.linalg.det(cov+np.eye(dim)*0.01) #協(xié)方差矩陣的行列式
    covinv = np.linalg.inv(cov+np.eye(dim)*0.01) #協(xié)方差矩陣的逆
    xdiff = x - mean
    #概率密度
    prob = 1.0/np.power(2*np.pi,1.0*dim/2)/np.sqrt(np.abs(covdet))*np.exp(-1.0/2*np.dot(np.dot(xdiff,covinv),xdiff))
    return prob


#獲取初始協(xié)方差矩陣
def getconvs(data,K):
    convs = [0]*K
    for i in range(K):
        # 初始的協(xié)方差矩陣源自于原始數(shù)據(jù)的協(xié)方差矩陣，且每個(gè)簇的初始協(xié)方差矩陣相同
        convs[i] = np.cov(data.T)  
    return convs


def isdistinct(means,criter=0.03): #檢測(cè)初始中心點(diǎn)是否靠得過(guò)近
    K = len(means)
    for i in range(K):
        for j in range(i+1,K):
            if criter > np.linalg.norm(means[i]-means[j]):
                return 0       
    return True
#獲取初始聚簇中心
def getmeans(data,K,criter):
    means = [0]*K
    dim  = np.shape(data)[1]
    minmax = [] #各個(gè)維度的極大極小值
    for i in range(dim):
        minmax.append(np.array([min(data[:,i]),max(data[:,i])]))

    while True:
        #生成初始點(diǎn)的坐標(biāo)
        for i in range(K):
            means[i] = []
            for j in range(dim):
                means[i].append(np.random.random()*(minmax[j][1]-minmax[j][0])+minmax[j][0])  
            means[i] = np.array(means[i])

        if isdistinct(means,criter):
            break  
    return means

# k-means算法的實(shí)現(xiàn)函數(shù)。
#用K-means算法輸出的聚類中心，作為高斯混合模型的輸入
def kmeans(data,K):
    N = np.shape(data)[0]#樣本數(shù)目
    dim = np.shape(data)[1] #維度

    means = getmeans(data,K,criter=15)
    means_old = [np.zeros(dim) for k in range(K)]


    while np.sum([np.linalg.norm(means_old[k]-means[k]) for k in range(K)]) > 0.01:

        means_old = cp.deepcopy(means)

        numlog = [0]*K
        sumlog = [np.zeros(dim) for k in range(K)]
        for n in range(N):
            distlog = [np.linalg.norm(data[n]-means[k]) for k in range(K)]
            toK = distlog.index(np.min(distlog))

            numlog[toK] += 1
            sumlog[toK] += data[n]

        for k in range(K):
            means[k] = 1.0/numlog[k]*sumlog[k]
    return means

#對(duì)程序結(jié)果進(jìn)行可視化，注意這里的K只能取2，否則該函數(shù)運(yùn)行出錯(cuò)
def visualresult(data,gammas,K):
    N = np.shape(data)[0]#樣本數(shù)目
    dim = np.shape(data)[1] #維度

    minmax = [] #各個(gè)維度的極大極小值
    xy = []
    n=200
    for i in range(dim):
        delta = 0.05*(np.max(data[:,i])-np.min(data[:,i]))
        xy.append(np.linspace(np.min(data[:,i])-delta,np.max(data[:,i])+delta,n))
    xx,yy = np.meshgrid(xy[0], xy[1])
    zz = np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            zz[i][j] = np.sum(gaussian(np.array([xx[i][j],yy[i][j]]),means[k],convs[k]) for k in range(K))
    gci = plt.imshow(zz,origin='lower',alpha = 0.8) # 選項(xiàng)origin='lower' 防止tuixan圖像顛倒
    plt.xticks([0,len(xy[0])-1],[xy[0][0],xy[0][-1]])
    plt.yticks([0,len(xy[1])-1],[xy[1][0],xy[1][-1]])

    for i in range(N):
        if gammas[i][0] >0.5:
            plt.plot((data[i][0]-np.min(data[:,0]))/(xy[0][1]-xy[0][0]),(data[i][1]-np.min(data[:,1]))/(xy[1][1]-xy[1][0]),'r.')
        else:
            plt.plot((data[i][0]-np.min(data[:,0]))/(xy[0][1]-xy[0][0]),(data[i][1]-np.min(data[:,1]))/(xy[1][1]-xy[1][0]),'k.')

    deltax = xy[0][1]-xy[0][0]
    deltay = xy[1][1]-xy[1][0]

    plt.plot((means[0][0]-xy[0][0])/deltax,(means[0][1]-xy[1][0])/deltay,'*r',markersize=15)
    plt.plot((means[1][0]-xy[0][0])/deltax,(means[1][1]-xy[1][0])/deltay,'*k',markersize=15)

    plt.title(u'高斯混合模型圖',{'fontname':'STFangsong','fontsize':18})
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
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

高斯混合模型的主函數(shù)

N = np.shape(data)[0]#樣本數(shù)目
dim = np.shape(data)[1] #維度
K = 2 # 聚簇的個(gè)數(shù)

means = kmeans(data,K)

convs = getconvs(data,K)

pis = [1.0/K]*K
gammas = [np.zeros(K) for i in range(N)] #*N 注意不能用 *N，否則N個(gè)array只指向一個(gè)地址

loglikelyhood = 0
oldloglikelyhood = 1

while np.abs(loglikelyhood - oldloglikelyhood)> 0.0001:
    oldloglikelyhood = loglikelyhood


    # E_step
    for n in range(N):
        respons = [pis[k]*gaussian(data[n],means[k],convs[k]) for k in range(K)]

        sumrespons = np.sum(respons)
        for k in range(K):
            gammas[n][k] = respons[k]/sumrespons

    # M_step
    for k in range(K):
        nk = np.sum([gammas[n][k] for n in range(N)])
        means[k] = 1.0/nk * np.sum([gammas[n][k]*data[n] for n in range(N)],axis=0)

        xdiffs = data - means[k]
        convs[k] = 1.0/nk * np.sum([gammas[n][k]*xdiffs[n].reshape(dim,1)*xdiffs[n] for n in range(N)],axis=0)
        pis[k] = 1.0*nk/N

    # 計(jì)算似然函數(shù)值
    loglikelyhood =np.sum( [np.log(np.sum([pis[k]*gaussian(data[n],means[k],convs[k]) for k in range(K)])) for n in range(N) ])
    #print means
    #print loglikelyhood
    #print '=='*10

visualresult(data,gammas,K)
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
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

4.高斯混合模型聚簇效果圖

5.參考文獻(xiàn)：

• K-means聚類和EM思想
  http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006910.html
• （EM算法）The EM Algorithm
  http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html
• 從決策樹學(xué)習(xí)談到貝葉斯分類算法、EM、HMM
  http://blog.csdn.net/v_july_v/article/details/7577684
• EM算法學(xué)習(xí)(Expectation Maximization Algorithm)
  http://www./XUHV3.clip
• EM算法——理論與應(yīng)用
  http://blog.sina.com.cn/s/blog_a8fead9b01014p6k.html
• EM算法 一個(gè)簡(jiǎn)單的例子
  http://blog.sina.com.cn/s/blog_a7da5cda010158b3.html
• 高斯混合模型（GMM）
  http://www.cnblogs.com/mindpuzzle/archive/2013/04/24/3036447.html
• 高斯混合模型
  http://www.cnblogs.com/zhangchaoyang/articles/2624882.html
• CS229 Lecture notes,Andrew Ng講義
  http://cs229./notes/cs229-notes7b.pdf
• https://github.com/lzhang10/maxent/blob/master/doc/manual.pdf
• http://nlp./software/classifier.shtml
• https://github.com/juandavm/em4gmm
  用最大期望算法求解高斯混合模型，dat文件夾里的壓縮文件是數(shù)據(jù)
• http:///python-unsupervised-learning-using-em-algorithm-implementation/
  Python實(shí)現(xiàn)的期望最大化算法，數(shù)據(jù)是鳶尾花數(shù)據(jù)