高斯混合模型

1、什麼是高斯混合模型（GMM）

　高斯混合模型（Gaussian Mixed Model）指的是多個高斯分佈函數的線性組合，一般用於解決同一集合下的數據包含多個不一樣的分佈的狀況，如解決分類狀況

　以下圖，明顯分紅兩個聚類。這兩個聚類中的點分別經過兩個不一樣的正態分佈隨機生成而來。若是隻用一個的二維高斯分佈來描述圖中的數據。這顯然不太合理，畢竟肉眼一看就應該分紅兩類

 

　但使用兩個二維高斯分佈來描述圖中的數據，將兩個二維高斯分佈\(N(\mu_1,\sum_1),N(\mu_2,\sum_2)\)作線性組合，用線性組合後的分佈來描述整個集合中的數據。這就是高斯混合模型（GMM）

2、GMM原理

　設有隨機變量\(X\)，則混合高斯模型能夠這樣表示： 

\(p(x) = \sum_{k=1}^{K}\pi_kN(x|\mu_k,\sum_k)\)

　其中\(N(x|\mu_k,\sum_k)\)稱爲混合模型中的第k個份量（component）。如前面圖中的例子，有兩個聚類，能夠用兩個二維高斯分佈來表示，那麼份量數K=2，\(\pi_k\)是混合係數（mixture coefficient），且知足：

\( \sum_{k=1}^{K}\pi_k,0\leq\pi_k\leq 1 \)，其中\(\pi_k\)至關於每一個高斯分佈的權重

 

　GMM用於聚類時，假設數據服從混合高斯分佈（Mixture Gaussian Distribution），例如上圖的例子，很明顯有兩個聚類，能夠定義K=2，那麼對應的GMM形式以下：

\( p(x) = \pi_1N(x|\mu_1,\sum_1) + \pi_2N(x|\mu_2,\sum_2) \)

 

　上式中未知參數有六個：\((\pi_1,\mu_1,\sum_1, \pi_2,\mu_2,\sum_2)\)，GMM聚類時分爲兩步，第一步是隨機地在這K個份量中選一個，每一個份量被選中的機率即爲混合係數\(\pi_k\)，能夠設定\(\pi_1=\pi_2=0.5\)，表示每一個份量被選中的機率是0.5，即從中抽出一個點，這個點屬於第一類的機率和第二類的機率各佔一半。但實際應用中事先指定\(\pi_k\)的值是很笨的作法：當從圖中的集合隨機選取一個點，怎麼知道這個點是來自\(N(x|\mu_1,\sum_1)\)仍是\(N(x|\mu_2,\sum_2)\)呢？

　換言之怎麼根據數據自動肯定\(\pi_1,\pi_2\)的值？這就是GMM參數估計的問題。要解決這個問題，可使用EM算法。經過EM算法，咱們能夠迭代計算出GMM中的參數：\(\pi_k,\mu_k,\sum_k\)

 

　咱們目的是找到這樣一組參數，它所肯定的機率分佈生成這些給定的數據點的機率最大，而這個機率實際上就等於：\(\prod_{i=1}^{N}p(x_i)\)，咱們把這個乘積稱做似然函數 (Likelihood Function)。咱們一般會對其取對數，把乘積變成加和：\(\sum_{i=1}^{N}logp(x_i)\)，獲得 log-likelihood function 。接下來將這個函數最大化（一般的作法是求導並令導數等於零，而後解方程），咱們就認爲這是最合適的參數，這樣就完成了參數估計的過程。下面讓咱們來看一看 GMM 的 log-likelihood function ：

\(\sum_{i=1}^{N}log\{\sum_{i=1}^{K}\pi_kN(x_i|\mu_k,\sum_k)\}\)

 

　因爲上式太複雜，無法直接用求導求得最大值。爲了解決這個問題，咱們採起以前從 GMM 中隨機選點的辦法：分紅兩步

 

　　1.估計每一個數據由每一個 Component 生成的機率（並非每一個 Component 被選中的機率）：對於每一個數據\(x_i\)來講，它由第\(k\)個 Component 生成的機率爲：

　　　因爲式子裏的\(\mu_k,\sum_k\)也是須要咱們估計的值，咱們採用迭代法，在計算\(\gamma(i,k)\)的時候咱們假定\(\mu_k,\sum_k\)均已知，咱們將取上一次迭代所得的值（或者初始值）

 

　　2.估計每一個Component 的參數：如今咱們假設上一步中獲得的\(\gamma(i,k)\)就是正確的"數據\(x_i\)由Component \(k\)生成的機率"，亦能夠當作該Component 在生成這個數據上所作的貢獻，或者說，咱們能夠看做\(x_i\)這個數據其中有\(\gamma(i,k)x_i\)這部分是由Component \(k\)所生成的，集中考慮全部的數據點，如今實際上能夠看做Component 生成了\( \gamma(1,k)x_1,...,\gamma(N,k)x_N \)這些點，因爲每一個 Component 都是一個標準的 Gaussian 分佈，能夠很容易求出最大似然所對應的參數值： 

　　　其中\(N_k = \sum_{i=1}^{N}\gamma(i,k)\)，而且\(\pi_k\)也瓜熟蒂落地能夠估計\(N_k/N\)

 

　　3.重複迭代前面兩步，直到似然函數的值收斂爲止

 

3、代碼實例

　參考：http://www.javashuo.com/article/p-sqhevuln-cx.html

'''
此示例程序隨機從4個高斯模型中生成500個2維數據，真實參數：
混合項w=[0.1，0.2，0.3，0.4]
均值u=[[5，35]，[30，40]，[20，20]，[45，15]]
協方差矩陣∑=[[30，0]，[0，30]]
而後以這些數據做爲觀測數據，根據EM算法來估計以上參數
此程序未估計協方差矩陣

'''

import math
import copy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


iter_num = 1000
N = 500
k = 4
probility = np.zeros(N)

u1 = [5, 35]
u2 = [30, 40]
u3 = [20, 20]
u4 = [45, 15]

sigma = np.matrix([[30, 0], [0,30]])
alpha = [0.1, 0.2, 0.3, 0.4]

#生成隨機數據，4個高斯模型
def generate_data(sigma, N, mu1, mu2, mu3, mu4, alpha):
          global X#可觀測數據集
          X = np.zeros((N, 2))#X：2維數據，N個樣本
          X = np.matrix(X)

          global mu#隨機初始化mu1，mu2，mu3，mu4
          mu = np.random.random((4, 2))
          mu = np.matrix(mu)

          global excep#指望：第i個樣本屬於第j個模型的機率的指望
          excep = np.zeros((N, 4))

          global alpha_#初始化混合項係數
          alpha_ = [0.25, 0.25, 0.25, 0.25]

          #np.random.multivariate_normal()：用於根據實際狀況生成一個多元正態分佈矩陣
          for i in range(N):
                    if np.random.random(1) < 0.1:
                              X[i, :] = np.random.multivariate_normal(mu1, sigma, 1)
                    elif 0.1 <= np.random.random(1) < 0.3:
                              X[i, :] = np.random.multivariate_normal(mu2, sigma, 1)
                    elif 0.3 <= np.random.random(1) < 0.6:
                              X[i, :] = np.random.multivariate_normal(mu3, sigma, 1)
                    else:
                              X[i, :] = np.random.multivariate_normal(mu4, sigma, 1)
          print('可觀測數據：\n', X)
          print('初始化的mu1, mu2, mu3, mu4：', mu)

def e_step(sigma, k, N):
          global X
          global mu
          global excep
          global alpha_

          for i in range(N):
                    denom = 0
                    for j in range(0, k):
                              denom += alpha_[j] * math.exp(-0.5 * (X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))#分母

                    for j in range(0, k):
                              numer = math.exp(-0.5 * (X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))#分子
                              excep[i, j] = alpha_[j] * numer / denom#求指望

          print('隱藏變量：\n', excep)
                    
def m_step(k, N):
          global excep
          global X
          global alpha_

          for j in range(0, k):
                    denom = 0#分母
                    numer = 0#分子
                    for i in range(N):
                              numer += excep[i, j] * X[i, :]
                              denom += excep[i, j]

                    mu[j, :] = numer / denom #求均值
                    alpha_[j] = denom / N#求混合項係數

generate_data(sigma, N, u1, u2, u3, u4, alpha)

#迭代計算
for i in range(iter_num):
          err = 0#均值偏差
          err_alpha = 0#混合係數偏差

          Old_mu = copy.deepcopy(mu)
          Old_alpha = copy.deepcopy(alpha_)

          e_step(sigma, k, N)#E 步
          m_step(k, N)#M 步

          print("迭代次數:", i+1)
          print("估計的均值:", mu)
          print("估計的混合項係數:", alpha_)
          for z in range(k):
                    err += (abs(Old_mu[z, 0] - mu[z, 0]) + abs(Old_mu[z, 1] - mu[z, 1]))#計算偏差
                    err_alpha += abs(Old_alpha[z] - alpha_[z])

          if (err <= 0.001) and (err_alpha < 0.001):#達到精度退出迭代
                    print(err, err_alpha)
                    break
          
#可視化結果，畫生成的原始數據
plt.subplot(221)
plt.scatter(X[:,0].tolist(), X[:,1].tolist(),c='b', s=25, alpha=0.4, marker='o')
plt.title('random generated data')

#畫分類好的數據
plt.subplot(222)
plt.title('classified data through EM')
order = np.zeros(N)
color = ['b', 'r', 'k', 'y']
for i in range(N):
          for j in range(k):
                    if excep[i, j] == max(excep[i, :]):
                              order[i] = j#選出X[i, :]屬於第幾個高斯模型
                    probility[i] += alpha_[int(order[i])] * math.exp(-0.5 * (X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/(np.sqrt(np.linalg.det(sigma))*2*np.pi)#計算混合高斯分佈

          plt.scatter(X[i, 0], X[i, 1], c=color[int(order[i])], s=25, alpha=0.4, marker='o')#繪製分類後的散點圖

#繪製三維圖像
ax = plt.subplot(223, projection='3d')
plt.title('3d view')
for i in range(N):
          ax.scatter(X[i,0], X[i, 1], probility[i], c=color[int(order[i])])
plt.show()