標籤（空格分隔）： 王小草機器學習筆記

EM演算法的英文全稱是Expectation Maximization Algorithm，也就是求期望最大化，也就是我們常說的目標函式求最大值的演算法。EM演算法，直觀的說，就是有一堆未知的資料（比如一些特徵值），這些資料可能來自於不同的類別，而你想知道的是每一個數據都來自於哪個類別，並且知道來自於這個類別的概率是多少。而在EM演算法看來，每一個類別中的資料必然服從了某個固有的分佈（如二項分佈，正態分佈等），只要尋找它密度函式的引數，就能知道資料的分佈，所以EM使用了極大似然估計的方法去估計了各個分佈的引數，從而進行了無監督地聚類。

本文首先會回顧和複習一下EM演算法中要涉及到的知識：Jensen不等式和最大似然估計；
然後介紹EM演算法，並且會用混合高斯GMM來作為例子用EM演算法求解。
最後，會給出一些實現EM演算法的python例項。

1. 回顧

1.1 回顧Jensen不等式

Jensen不等式這個名字有點陌生，但是如果眼睜睜地看到這個不等式，你肯定會覺得特別特別熟悉，並且鼻尖傳來陣陣數學考試卷的味道，它會出現在選擇題填空題或者後面的大分推導題中。萬萬沒想到，今天還會在那麼重要的演算法中遇到。

假設有兩個變數x1,x2,有一個函式f(),並且函式f是凸函式，那麼就肯定有以下不等式成立：

f(θx1 + (1-θ)x2) <= θf(x1) + (1-θ)f(x2)

在二維座標中表示如下

同理，若有多個x，多個θ，只要滿足f是凸函式，並且
θ1，θ2，…θk >= 0;且θ1+θ2+…+θk = 1
那麼下面不等式肯定成立：

上面講的是離散型的變數，針對連續型的變數，Jesen不等式也是成立的。

中θ是隨機變數，分別與xi相乘後相加，其實就是在求x的期望值，那麼以上任何形式的不等式，都可以表示成如下：

1.2 回顧K-means演算法

之所以要在EM演算法中回顧K-means，是因為在迭代k-means中其實就是不斷取均值點求期望的過程，於EM迭代求最大化期望類似。

先初始選擇k個簇中心，然後根據各個樣本點到中心的距離來分割，再對分割後的類找到這個類中所有樣本點的中心點（均值點）作為新的簇中心，然後再進行聚類，再找新的中心，如此迴圈，直到滿足終止條件而停止。

這裡，通過各個樣本的均值來確定新的中心點其實就是在求期望，其過程的大致思路其實與我們接下去要講的EM演算法是類似的。

k-means能夠非常方便得將樣本分成若干簇，但是由一個巨大的問題是，我們無法知道這個樣本點屬於這個簇的概率，我們只知道屬於或不屬於的布林值，如此一來，就沒那麼好玩了。

那麼如果想知道概率的話，就得去找到與樣本分佈最接近的概率分佈模型，要得到概率分佈模型，就得先將模型中的引數估計出來。EM演算法所實現的就是這個功能。

2.最大似然估計求引數

2.1 小引子

先來看一個簡單的例子。
簡單的例子一般都離不開拋硬幣來。

假設我們現在不知道拋硬幣出現正反面的概率，然後設每拋一次出現正面的概率是p，這個p是我們想求得的。
現在我做一個實驗，將硬幣拋是10次，然後記錄結果為：正正反正正正反反正正。
最大似然估計就是去求出現以上這10次結果的概率最大時候，去估計p的概率。
因為每次拋硬幣是一個獨立事件，所以每一次拋硬幣的概率是可以相乘的，於是以上10次結果發生的概率可以寫成：

要求的P最大時p的值，其實就是對以上等式求目標函式最優化的過程，最後可以求出p=0.7。

當然，你肯定會義正言辭地說，不對！拋硬幣的概率誰不知道呢，正反兩面出現的概率都是0.5呀！
這是因為我們這個實驗中只拋了10次，樣本量太小存在的誤差會偏大，如果拋100次，1000次，10000次，樣本中正反兩面出現的次數會越來越趨於1:1，那麼求出來的p值肯定也會更接近與0.5了。

當然，你肯定還會義正言辭得說，可是這個極大似然估計有什麼用嗎，拋硬幣的概率我本來就知道。嗯，對，拋硬幣只是一個例子。假如我收集了上海9月份30天的天氣資料，想知道9月份的上海下雨的概率p有多大；再假如我記錄了今年我從東昌路上2號線有沒有座位的資料，想知道上車後有座位的概率p是多少；再假如我有所有進入購物網站的行為資料樣本，想知道首次進來的人會購買下單的概率；再假如二號店若推出新品促銷的廣告，使用者看到會點選進入的概率…等等。這些概率我想你應該不知道，但作為商家也許會非常渴望知道。

2.2 二項分佈的最大似然估計

拋硬幣其實是一個二項分佈，它有兩個值，概率分別為p和1-p。
假設投擲N次硬幣，出現正面朝上次數是n，反面朝上的次數是N-n。
並且設正面朝上的概率是p。
現在使用似然函式來求目標函式的最優化，為了計算方便，我們將函式取對數來求最大值，成為“對數似然函式”。
目標函式如下：

對目標函式中的p求導數，最後得到p = n/N

以上就是用最大似然函式估計二項分佈引數的過程。

2.3 高斯分佈的最大似然估計

現在我們來看一看高斯分佈。
若給定一組樣本X1,X2,X3…Xn，已知他們是來自於高斯分佈N（μ，σ)，即符合均值為μ，標準差為σ的正太分佈。要根據這些樣本點的分佈去估計這個正態分佈的引數μ，σ。

1.首先，要知道高斯分佈的概率密度函式：

裡面有兩個引數，分別就是μ，σ。也是我們要求的引數。

2.要得到樣本點那樣分佈的概率，假設每個樣本都是獨立的，所有總的概率就是每個樣本的概率的乘積：

3.對上面的等式進行對數似然函式的化簡：

4.得到化簡後的目標函式：

現在就是對這個目標函式求最大值時的引數μ，σ的值

5.將目標函式分別對引數μ，σ進行求導，就能求出μ，σ的公式：

以上就是用最大似然函式估計高斯分佈引數的過程。

3.EM演算法

經過了冗長的鋪墊終於到了本文的重點了–EM演算法。在講EM演算法乏味難懂的定理前，我們先用高斯混合模型來走一遍它的引數估計

3.1 直觀理解猜測GMM的引數估計

已知一個學校的所有學生的身高樣本（X1,X2,X3…Xn），並且男生和女生都分別服從N（μ男，σ男）和N（μ女，σ女）的高斯分佈。目的是要求出μ男，σ男，μ女，σ女這四個引數。

來來來理一下這個題目，與上文中的例子不同了，現在同時有兩個高斯分佈混合在一起，我們要去求出兩個高斯分佈各自的均值與標準差引數。那麼問題來了…我怎麼知道某個樣本屬於男高斯還是女高斯啊，現在我們只知道所有的身高值，並不知道這些身高背後的性別呀。恩恩，這是一個混合高斯模型，簡稱GMM（隨機變數是由2個高斯分佈混合而成）。在GMM中，有一個隱藏的隨機變數我們沒法看到，那就是性別，因為無法知道性別的概率，也就無法知道某個樣本屬於哪個高斯分佈的了。所以這個隱藏的概率π男，π女也是我們需要估計的，它表示某個樣本屬於某個高斯分佈的概率。

將上面的例子擴充套件地表示出來：隨機變數X是由K個高斯分佈混合而成，取各個高斯分佈的概率為π1，π2，π3…πk；第i個高斯分佈的均值為μi,方差為Σi。若觀測到隨機變數X的一系列樣本X1,X2,X3…Xn，試估計引數π，μ，Σ。

1.建立目標函式
同樣的，我們使用對數似然函式來建立目標函式

由於對數函式裡面又有加和，無法直接用求導解方程的方法來求最大值。為了解決這個問題，接下來分兩步走。

2.第一步：估計資料來自哪個組（也就是說來自哪個高斯分佈）
首先我們根據經驗隨便給定π，μ，Σ的先驗值。
然後求r(i,k),表示，樣本i 由高斯分佈k生成的概率。公式如下：

根據這個公式，我們可以求得樣本X1分別屬於K1,K2,K3..的概率

3.第二步：估計每個高斯分佈中的引數
對於高斯分佈k來說，它說生成的點可看成是{r(i,k)*xi|i-=,1,2,3..N}，就是所有原來的樣本點乘以它屬於高斯分佈k的概率，從而得到了新的樣本點，這些樣本點應是服從高斯分佈k的。因此高斯分佈k的樣本個數不再是原來的N，而是所有樣本屬於它的概率的加和：
同理，π，μ，Σ都可以因此重新計算：

4.重複第一步，第二步
我們用先驗的π，μ，Σ計算出來新的π，μ，Σ，於是又可以利用新的π，μ，Σ去計算新的r(i,k)，再得出又一輪新的π，μ，Σ，如此迴圈往復，這些引數會慢慢收斂，直到前後兩次的差值小於預先設定的閥值，就停止迭代，最後一次算出的π，μ，Σ就可以當成最終的估計值啦。

3.2 EM演算法的提出

問題的提出：
假設有訓練樣本X1,X2,X3…Xm，包含m個獨立樣本，希望從中找出p(x,z)的引數。

通過似然估計建立目標函式：

x是樣本點，z是隱藏引數（就是上文中的π），θ是顯引數（就是上文高斯分佈中的μ，σ）

z是隱藏隨機變數，不方便直接找到引數估計。所以使用策略：計算l(θ)的下界，求出該下界的最大值，重複該過程知道收斂，下圖中綠線所經過的點是下界與l(θ)相等的點，求下界的最大值，往往會更接近目標函式的最大值。

設Qi是z的某一個分佈，那麼目標函式可以轉換：

log函式是一個凹函式，要求的一個點與它相等，只能使得是一個常數，也就是：

進一步分析：

到後來，居然發現Q(z)其實就是z關於樣本xi，引數θ的條件概率

由此，可以歸納出EM演算法的整體框架：

3.3 理論公式推導GMM

還是3.1中的例子：
已知一個學校的所有學生的身高樣本（X1,X2,X3…Xn），並且男生和女生都分別服從N（μ男，σ男）和N（μ女，σ女）的高斯分佈。目的是要求出μ男，σ男，μ女，σ女這四個引數。

E-STEP:

M-STEP:
將多項式和高斯分佈等引數帶入

對均值求偏導，求均值：

上式等於0解均值為：

對方差求偏導，求方差：

對多項式中的引數：
考察m-step中的目標函式，對於φ，刪除常數項

得到：

由於多項式的概率和為1，建立拉格朗日方程：

求偏導

於是我們仍然可以得到引數的估計方程：

與3.1相比，得到了相同的結果！

5. Python實現

5.1 EM演算法估算GMM均值

這個例子是我網上找的，註釋不多，所以我將每句程式碼都寫了註釋以方便理解。

第一個方法ini_data是生成了一組由兩個高斯分佈產生的樣本資料，他們的方差相同，均值不同。
第二個方法e_step是EM演算法的第1步，先初始化引數，然後根據已知的引數去估計概率π
第三個方法M-step是EM演算法的第2步，根據e-step中計算出來的π去估計新的μ和σ
第四個方法run是去不斷迭代迴圈e-step和m-step直至收斂求得最終的三個引數的估計值

#!/usr/bin/python
# -*- coding:utf-8 -*-

import math
import copy
import numpy as np
import matplotlib.pyplot as plt

is_debug = False


# # 建立兩個正態分佈的資料，方差已知並同，均值不同，k是高斯分佈的個數，n是樣本數
def ini_data(sigma, mu1, mu2, k, n):
    global x
    global mu
    global expectations
    # 建立1*n維的0向量
    x = np.zeros((1, n))
    # 隨機生成兩個0-1間的隨機數
    mu = np.random.random(2)
    # 生成n*k維的0向量
    expectations = np.zeros((n, k))
    # 遍歷
    for i in xrange(0, n):
        # 如果生成一個隨機數大於0.5
        if np.random.random(1) > 0.5:
            # 那麼x中第0行第i列由第1個高斯分佈生成
            x[0,i] = np.random.normal()*sigma + mu1
        # 如果小於0.5
        else:
            # 則由第2個高斯分佈生成隨機數
            x[0,i] = np.random.normal()*sigma + mu2
    if is_debug:
        print "***"
        print u"初始觀測值資料x:"
        print x


# # 計算E[Zij]隱藏引數
def e_step(sigma, k, n):
    global expectations
    global mu
    global x
    # 對每個樣本做遍歷
    for i in xrange(0, n):
        # 初始化新的樣本數
        denom = 0
        # 對該樣本在每個分佈中做遍歷
        for j in xrange(0, k):
            # 累加計算所有樣本點在該分佈中的數
            denom += math.exp((-1/(2 * float(sigma ** 2))))*(float(x[0,i] - mu[j]) ** 2)
        for j in xrange(0, k):
            # 計算該樣本點在該分佈中的數
            numer = math.exp((-1/(2 * float(sigma ** 2))))*(float(x[0,i] - mu[j]) ** 2)
            # 計算該樣本點來自該分佈的概率r(i,k)
            expectations[i,j] = numer / denom
    if is_debug:
        print "***"
        print u"隱藏變數E(Z)"
        print expectations


# # M-step:求最大化E(Zij)時的引數mu
def m_step(k, N):
    global expectations
    global x
    # 遍歷每個分佈
    for j in xrange(0, k):
        # 初始化新的樣本點的值
        numer = 0
        # 初始化新的樣本數量
        denom = 0
        # 對每個新樣本做遍歷
        for i in xrange(0, N):
            # 累加所有樣本的值
            numer += expectations[i, j] * x[0, i]
            # 累加所有樣本的數量
            denom += expectations[i, j]
        # 計算新的均值
        mu[j] = numer / denom


# # 迭代收斂
def run(sigma, mu1, mu2, k, n, iter_num, epsilon):
    # 生成由兩個正態分佈組成的資料
    ini_data(sigma, mu1, mu2, k, n)
    print u"初始化均值：", mu
    # 根據預先設定的迭代次數迴圈
    for i in range(iter_num):
        # 將現有的均值賦值給老的均值變數
        old_mu = copy.deepcopy(mu)
        # 計算隱變數的期望
        e_step(sigma, k, n)
        # 求目標函式最大值時的均值
        m_step(k, n)
        print i, mu
        # 如果新舊均值的差值小於設定的閥值則終止迴圈
        if sum(abs(mu - old_mu)) < epsilon:
            break


if __name__ == '__main__':
    run(6, 40, 20, 2, 10000, 10, 0.01)
    plt.hist(x[0, :], 50)
    plt.show()

最後一個方法將這組資料在二維座標上畫出來了，影象如下，確實可以很清晰地看出有兩個正態分佈，他們的均值不同，但方差類似。

打印出迭代了10次的均值結果：

初始化的均值是隨機的，非常小，與我們實際的均值差別非常大，在經過第一次迭代後，就已經有了明顯的改善了，接下來每次迭代的結果都越來越接近真實的均值40,20了，可以看到如果迭代次數增加到100，或1000,EM估計出來的結果應該會非常接近40和20的。

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from scipy.stats import multivariate_normal
import matplotlib as mpl
import matplotlib.pyplot as plt


mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False


# # 生成來自兩個高斯分佈的資料
if __name__ == '__main__':
    # 第一個分佈的均值
    mu1 = (0, 0, 0)
    # 建立一個維度為3的單位矩陣，作為方差
    cov1 = np.identity(3)
    # 根據均值與方差隨機生成多元的正太分佈的資料
    data1 = np.random.multivariate_normal(mu1, cov1, 100)
    # 第二個分佈的均值
    mu2 = (2, 2, 1)
    # 設定方差
    cov2 = np.identity(3)
    # 根據均值與方差隨機生成第二個多元的正太分佈的資料
    data2= np.random.multivariate_normal(mu2, cov2, 100)
    # 將兩組資料union
    data = np.vstack((data1, data2))

    # 設定迭代次數
    num_iter = 100
    # 樣本資料的大小，n是樣本數的個數，d是維度，這裡是3維
    n, d = data.shape
    # 初始化引數，隨機設定即可
    mu1 = np.random.standard_normal(d)  # 從三維的標準正態分佈中隨機取一個點作為初始均值
    mu2 = np.random.standard_normal(d)
    sigma1 = np.identity(d)  # 建立d維的單位矩陣作為初始化的方差
    sigma2 = np.identity(d)
    pi = 0.5 # 再拍腦門決定一下隱引數的值

    # # 進行EM演算法
    for i in range(num_iter):
        # E-step
        # 根據均值方差建立兩個正太分佈
        norm1 = multivariate_normal(mu1, cov1)
        norm2 = multivariate_normal(mu2, cov2)
        # 分佈計算每個樣本點由兩個分佈產生的概率
        tau1 = pi*norm1.pdf(data)
        tau2 = (1 - pi)*norm2.pdf(data)
        # 計算由第一個分佈產生的概率
        gamma = tau1/(tau1 + tau2)

        # M-step
        mu1 = np.dot(gamma, data)/sum(gamma)
        mu2 = np.dot((1-gamma), data)/sum(1-gamma)
        sigma1 = np.dot(gamma * (data - mu1).T, (data - mu1))/sum(gamma)
        sigma2 = np.dot((1-gamma) * (data - mu2).T, (data - mu2))/sum(1-gamma)
        pi = sum(gamma)/n
        # if i % 2 == 0:
        #     print i, ":\t",mu1, mu2

    print '類別概率：\t', pi
    print '均值：\t', mu1, mu2
    print '方差：\t', sigma1, sigma2

    g = GMM(n_components=2, covariance_type="full", n_iter=100)
    g.fit(data)
    print '類別概率:\t', g.weights_[0]
    print '均值:\n', g.means_, '\n'
    print '方差:\n', g.covars_, '\n'

    # 預測分類
    norm1 = multivariate_normal(mu1, sigma1)
    norm2 = multivariate_normal(mu2, sigma2)
    tau1 = norm1.pdf(data)
    tau2 = norm2.pdf(data)

    # figsize設定畫布的大小,facecolor設定畫布背景色為白色
    fig = plt.figure(figsize=(14, 7), facecolor='w')
    # 121表示畫1行2列中的第1個圖：原始資料的分佈圖
    ax = fig.add_subplot(121, projection='3d')
    # x,y,z座標是原始資料的第1,2,3列資料，c是顏色為藍色，s是三點的大小，marker是三點的形狀為圓，depthshade為顏色需不需要分深淺度
    ax.scatter(data[:, 0], data[:, 1], data[:, 2], c='b', s=30, marker='o', depthshade=True)
    # 設定座標的標記
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    # 設定影象的名稱與此題大小
    ax.set_title(u'origin data', fontsize=18)

    # 畫1行2列中的第2個圖：預測資料的分類圖
    ax = fig.add_subplot(122, projection='3d')
    # 提取出概率在分佈1中比分佈2中大的樣本點
    c1 = tau1 > tau2
    # 將這點畫在第2張圖中，顏色為紅色，形狀是小圓
    ax.scatter(data[c1, 0], data[c1, 1], data[c1, 2], c='r', s=30, marker='o', depthshade=True)
    # 同理，拿出在分佈2中概率大的樣本點
    c2 = tau1 < tau2
    # 用綠色的小三角表示
    ax.scatter(data[c2, 0], data[c2, 1], data[c2, 2], c='g', s=30, marker='^', depthshade=True)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title(u'EM演算法分類', fontsize=18)

    # 自動調整影象大小，使兩個圖形比較緊湊
    plt.tight_layout()
    # 畫出來吧，少年
    plt.show()

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from scipy.stats import multivariate_normal
from sklearn.cross_validation import train_test_split
from sklearn.mixture import GMM
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors


mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False


def expand(a, b):
    d = (b - a) * 0.05
    return a-d, b+d


if __name__ == '__main__':
    # 讀入資料（性別，身高，體重）
    data = np.loadtxt('D:\Documents\cc\python\data\HeightWeight.csv',
                      dtype = np.float, delimiter=',', skiprows=1)
    # 將資料分成兩部分y:性別標籤列，x:身高與體重2列,[1, ]表示第一部分取[0,1)列，第二部分取剩下的列
    y, x = np.split(data, [1, ], axis=1)
    # 將資料分成訓練集與測試集，random_state保證了隨機池相同
    x, x_test, y, y_test = train_test_split(x, y, train_size=0.6, random_state=0)
    # 用生成一個GMM，由兩個分佈組成，tol是閥值
    gmm = GMM(n_components=2, covariance_type='full', tol=0.0001, n_iter=100, random_state=0)
    # 分別取出身高與體重的最大值與最小值
    x_min = np.min(x, axis=0)
    x_max = np.max(x, axis=0)
    # 使身高與體重的訓練資料去擬合剛剛建立的gmm
    gmm.fit(x)
    # 列印均值，方差看一看
    print '均值 = \n', gmm.means_
    print '方差 = \n', gmm.covars_
    # 對訓練資料與測試資料都帶入gmm做預測
    y_hat = gmm.predict(x)
    y_test_hat = gmm.predict(x_test)
    # 以為gmm做預測是不會分兩個分佈的先後順序的，所以要保證是順序顛倒的時候女性仍然表示為0
    change = (gmm.means_[0][0] > gmm.means_[1][0])
    if change:
        z = y_hat == 0
        y_hat[z] = 1
        y_hat[~z] = 0
        z = y_test_hat == 0
        y_test_hat[z] = 0
        y_test_hat[~z] = 0
    # 求訓練與測試集的準確率
    acc = np.mean(y_hat.ravel() == y.ravel())
    acc_test = np.mean(y_test_hat.ravel() == y_test.ravel())
    acc_str = u'訓練集準確率：%.2f%%' % (acc * 100)
    acc_test_str = u'測試集準確率：%.2f%%' % (acc_test * 100)
    print acc_str
    print acc_test_str

    # 建立淺色的顏色（用來做背景）和深色的顏色（用來表示點）
    cm_light = mpl.colors.ListedColormap(['#FF8080', '#77E0A0'])
    cm_dark = mpl.colors.ListedColormap(['r', 'g'])
    # 分別取兩個特徵的最大值與最小值
    x1_min, x1_max = x[:, 0].min(), x[:, 0].max()
    x2_min, x2_max = x[:, 1].min(), x[:, 1].max()
    x1_min, x1_max = expand(x1_min, x1_max)
    x2_min, x2_max = expand(x2_min, x2_max)
    # 在畫布上畫上橫縱座標網格各500
    x1, x2 = np.mgrid[x1_min:x1_max:500j, x2_min:x2_max:500j]
    # 將兩個特徵的數鋪平成一列
    grid_test = np.stack((x1.flat, x2.flat), axis=1)
    # 用gmm做預測y值
    grid_hat = gmm.predict(grid_test)
    # 將預測出來的y值調整大小成與x1一致
    grid_hat = grid_hat.reshape(x1.shape)
    if change:
        z = grid_hat == 0
        grid_hat[z] = 1
        grid_hat[~z] = 0
    # 畫一個9*7的畫布，背景色為白色
    plt.figure(figsize=(9, 7), facecolor='w')
    # 用之前設定的淺顏色來做兩個分佈的區分割槽域
    plt.pcolormesh(x1, x2, grid_hat,cmap=cm_light)
    # 畫訓練樣本點，用圓表示，並使用之前設定的深色表示
    plt.scatter(x[:, 0], x[:, 1], s=50, c=y, marker='o', cmap=cm_dark, edgecolors='k')
    # 畫訓練樣本點，用三角形表示，並使用之前設定的深色表示
    plt.scatter(x_test[:, 0], x_test[:, 1], s=60, c=y_test, marker='^', cmap=cm_dark, edgecolors='k')
    # 奔跑吧，少年
    plt.show()

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from sklearn.mixture import GMM, DPGMM
from sklearn.cross_validation import train_test_split
import matplotlib as mpl
import matplotlib.colors
import matplotlib.pyplot as plt

mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False


def expand(a, b, rate = 0.05):
    d = (b - a) * rate
    return a-d, b+d


def accuracy_rate(y1, y2):
    acc1 = np.mean(y1 == y2)
    acc = acc1 if acc1 > 0.5 else 1-acc1
    return acc


if __name__ == '__main__':
    # 設定隨機池
    np.random.seed(0)
    # 設定協方差
    cov1 = np.diag((1, 2))
    # 設定生成隨機數的個數
    N1 = 500
    N2 = 300
    N = N1 + N2
    # 生成正態分佈的隨機數
    x1 = np.random.multivariate_normal(mean=(3, 2), cov=cov1, size=N1)
    # 建立一個數組
    m = np.array(((1, 1), (1, 3)))
    # 將x1與新陣列做點乘，目的是做一個轉換
    x1 = x1.dot(m)
    x2 = np.random.multivariate_normal(mean=(-1, 10), cov=cov1, size=N2)
    # union兩組隨機數
    x = np.vstack((x1, x2))
    # 設定它們的y值
    y = np.array([0]*N1 + [1]*N2)

    # 型別
    types = ('spherical', 'diag', 'tied', 'full')
    err = np.empty(4)
    bic = np.empty(4)
    # 對型別進行列舉，i是索引
    for i, type in enumerate(types):
        gmm = GMM(n_components=2, covariance_type=type, random_state=0)
        gmm.fit(x)
        err[i] = 1 - accuracy_rate(gmm.predict(x), y)
        bic[i] = gmm.bic(x)
    print '錯誤率：', err.ravel()
    print 'BIC：', bic.ravel()

錯誤率： [ 0.44625  0.3625   0.30625  0.     ]
BIC： [ 8291.31989604  8189.04587246  8290.20374814  6850.90774844]
均值 = 
[[ -0.98545235  10.07568825]
 [  4.88245339   8.69755024]]
方差 = 
[[[  0.89173153  -0.02570634]
  [ -0.02570634   1.95207306]]

 [[  2.86753624   6.6289323 ]
  [  6.6289323   17.97477745]]]