本文翻譯自https://jakevdp.github.io/PythonDataScienceHandbook/05.12-gaussian-mixtures.html

上一節(jié)中探討的k-means聚類模型簡單易懂，但其簡單性導(dǎo)致其應(yīng)用中存在實際挑戰(zhàn)条辟。具體而言黔夭，k-means的非概率特性及簡單地計算點與類蔟中心的歐式距離來判定歸屬，會導(dǎo)致其在許多真實的場景中性能較差羽嫡。本節(jié)本姥，我們將探討高斯混合模型(GMMs)，其可以看成k-means的延伸厂僧，更可以看成一個強有力的估計工具扣草，而不僅僅是聚類。

我們將以一個標(biāo)準(zhǔn)的import開始

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np

GMM的動機：k-means的弱點

我們看下k-means的缺陷颜屠，思考下如何提高聚類模型辰妙。正如上一節(jié)所示，給定簡單甫窟，易于分類的數(shù)據(jù)密浑，k-means能找到合適的聚類結(jié)果。
舉例而言粗井，假設(shè)我們有些簡單的數(shù)據(jù)點尔破，k-means算法能以某種方式很快地將它們聚類，跟我們?nèi)庋鄯直娴慕Y(jié)果很接近：

# Generate some data
from sklearn.datasets.samples_generator import make_blobs
X, y_true = make_blobs(n_samples=400, centers=4,
                       cluster_std=0.60, random_state=0)
X = X[:, ::-1] # flip axes for better plotting

# Plot the data with K Means Labels
from sklearn.cluster import KMeans
kmeans = KMeans(4, random_state=0)
labels = kmeans.fit(X).predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');

從直觀的角度來看浇衬，我可能期望聚類分配時懒构，某些點比其他的更確定：舉例而言，中間兩個聚類之間似乎存在非常輕微的重疊耘擂，這樣我們可能對這些數(shù)據(jù)點的分配沒有完全的信心胆剧。不幸的是，k-means模型沒有聚類分配的概率或不確定性的內(nèi)在度量（盡管可能使用bootstrap 的方式來估計這種不確定性）醉冤。為此秩霍，我們必須考慮泛化這種模型。
k-means模型的一種理解思路是蚁阳，它在每個類蔟的中心放置了一個圈（或者铃绒，更高維度超球面)，其半徑由聚類中最遠的點確定螺捐。該半徑充當(dāng)訓(xùn)練集中聚類分配的一個硬截斷：任何圈外的數(shù)據(jù)點不被視為該類的成員颠悬。我們可以使用以下函數(shù)可視化這個聚類模型：

from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist

def plot_kmeans(kmeans, X, n_clusters=4, rseed=0, ax=None):
    labels = kmeans.fit_predict(X)

    # plot the input data
    ax = ax or plt.gca()
    ax.axis('equal')
    ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)

    # plot the representation of the KMeans model
    centers = kmeans.cluster_centers_
    radii = [cdist(X[labels == i], [center]).max()
             for i, center in enumerate(centers)]
    for c, r in zip(centers, radii):
        ax.add_patch(plt.Circle(c, r, fc='#CCCCCC', lw=3, alpha=0.5, zorder=1))

kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X)

觀察k-means的一個重要發(fā)現(xiàn)矮燎，這些聚類模式必須是圓形的。k-means沒有內(nèi)置的方法來計算橢圓形或橢圓形的簇椿疗。因此漏峰，舉例而言糠悼，假設(shè)我們將相同的數(shù)據(jù)點作變換届榄，這種聚類分配方式最終變得混亂：

rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))

kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X_stretched)

E-M的推廣：高斯混合模型

高斯混合模型（GMM）試圖找到一個多維高斯概率分布的混合，以模擬任何輸入數(shù)據(jù)集器一。在最簡單的情況下课锌，GMM可用于以與k-means相同的方式聚類。

from sklearn.mixture import GMM
gmm = GMM(n_components=4).fit(X)
labels = gmm.predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');

但因為GMM包含概率模型盹舞，因此可以找到聚類分配的概率方式 - 在Scikit-Learn中产镐，通過調(diào)用predict_proba方法實現(xiàn)。它將返回一個大小為[n_samples, n_clusters]的矩陣踢步，用于衡量每個點屬于給定類別的概率：

probs = gmm.predict_proba(X)
print(probs[:5].round(3))

[[ 0.     0.     0.475  0.525]
 [ 0.     1.     0.     0.   ]
 [ 0.     1.     0.     0.   ]
 [ 0.     0.     0.     1.   ]
 [ 0.     1.     0.     0.   ]]

我們可以可視化這種不確定性癣亚，比如每個點的大小與預(yù)測的確定性成比例；如下圖获印，我們可以看到正是群集之間邊界處的點反映了群集分配的不確定性：

size = 50 * probs.max(1) ** 2  # square emphasizes differences
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size);

本質(zhì)上說述雾，高斯混合模型與k-means非常相似：它使用期望-最大化的方式，定性地執(zhí)行以下操作：

1. 選擇位置和形狀的初始猜想
2. 重復(fù)直到收斂
   1. E步驟：對于每個點，計算其屬于每個類別的概率權(quán)重
   2. M步驟：對于每個類別玻孟，使用E步算出的權(quán)重唆缴，根據(jù)所有數(shù)據(jù)點，更新其位置，規(guī)范化和形狀
      結(jié)果是，每個類別不是被硬邊界的球體界定近速，而是平滑的高斯模型悍手。正如在k-means的期望-最大方法一樣，這種算法有時可能會錯過全局最優(yōu)解，因此在實踐中使用多個隨機初始化。
      讓我們創(chuàng)建一個函數(shù)，通過基于GMM輸出霎匈，繪制橢圓來幫助我們可視化GMM聚類的位置和形狀：

from matplotlib.patches import Ellipse

def draw_ellipse(position, covariance, ax=None, **kwargs):
    """Draw an ellipse with a given position and covariance"""
    ax = ax or plt.gca()
    
    # Convert covariance to principal axes
    if covariance.shape == (2, 2):
        U, s, Vt = np.linalg.svd(covariance)
        angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
        width, height = 2 * np.sqrt(s)
    else:
        angle = 0
        width, height = 2 * np.sqrt(covariance)
    
    # Draw the Ellipse
    for nsig in range(1, 4):
        ax.add_patch(Ellipse(position, nsig * width, nsig * height,
                             angle, **kwargs))
        
def plot_gmm(gmm, X, label=True, ax=None):
    ax = ax or plt.gca()
    labels = gmm.fit(X).predict(X)
    if label:
        ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
    else:
        ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)
    ax.axis('equal')
    
    w_factor = 0.2 / gmm.weights_.max()
    for pos, covar, w in zip(gmm.means_, gmm.covars_, gmm.weights_):
        draw_ellipse(pos, covar, alpha=w * w_factor)

有了這個，我們可以看看四成分的GMM為我們的初始數(shù)據(jù)提供了什么：

gmm = GMM(n_components=4, random_state=42)
plot_gmm(gmm, X)

同樣送爸，我們可以使用GMM方法來擬合我們的拉伸數(shù)據(jù)集铛嘱；允許full的協(xié)方差，該模型甚至可以適應(yīng)非常橢圓形袭厂，伸展的聚類模式：

gmm = GMM(n_components=4, covariance_type='full', random_state=42)
plot_gmm(gmm, X_stretched)

這清楚地表明GMM解決了以前遇到的k-means的兩個主要實際問題墨吓。

選擇協(xié)方差類型

如果看了之前擬合的細節(jié)，你將看到covariance_type選項在每個中都設(shè)置不同嵌器。該超參數(shù)控制每個類簇的形狀的自由度肛真；對于任意給定的問題，必須仔細設(shè)置爽航。默認值為covariance_type =“diag”蚓让，這意味著可以獨立設(shè)置沿每個維度的類蔟大小，并將得到的橢圓約束為與軸對齊讥珍。一個稍微簡單和快速的模型是covariance_type =“spherical”历极，它約束了類簇的形狀，使得所有維度都相等衷佃。盡管它并不完全等效趟卸，其產(chǎn)生的聚類將具有與k均值相似的特征。更復(fù)雜且計算量更大的模型（特別是隨著維數(shù)的增長）是使用covariance_type =“full”氏义，這允許將每個簇建模為具有任意方向的橢圓锄列。
對于一個類蔟，下圖我們可以看到這三個選項的可視化表示：

GMM 作為密度估計

盡管GMM通常被歸類為聚類算法惯悠，但從根本上說它是一種密度估算算法邻邮。也就是說，GMM適合某些數(shù)據(jù)的結(jié)果在技術(shù)上不是聚類模型克婶，而是描述數(shù)據(jù)分布的生成概率模型筒严。
例如丹泉，考慮一下Scikit-Learn的make_moons函數(shù)生成的一些數(shù)據(jù)：

from sklearn.datasets import make_moons
Xmoon, ymoon = make_moons(200, noise=.05, random_state=0)
plt.scatter(Xmoon[:, 0], Xmoon[:, 1]);

如果我們嘗試用視為聚類模型的雙成分的GMM模擬數(shù)據(jù)，則結(jié)果不是特別有用：

gmm2 = GMM(n_components=2, covariance_type='full', random_state=0)
plot_gmm(gmm2, Xmoon)

但是如果我們使用更多成分的GMM模型鸭蛙，并忽視聚類的類別摹恨，我們會發(fā)現(xiàn)更接近輸入數(shù)據(jù)的擬合：

gmm16 = GMM(n_components=16, covariance_type='full', random_state=0)
plot_gmm(gmm16, Xmoon, label=False)

這里，16個高斯分布的混合不是為了找到分離的數(shù)據(jù)簇娶视，而是為了對輸入數(shù)據(jù)的整體分布進行建模晒哄。這是分布的一個生成模型，這意味著GMM為我們提供了生成與我們的輸入類似分布的新隨機數(shù)據(jù)的方法歇万。例如揩晴，以下是從這個16分量GMM擬合到我們原始數(shù)據(jù)的400個新點：

Xnew = gmm16.sample(400, random_state=42)
plt.scatter(Xnew[:, 0], Xnew[:, 1]);

GMM非常方便勋陪，可以靈活地建模任意多維數(shù)據(jù)分布贪磺。

多少components？

GMM是一種生成模型這一事實為我們提供了一種確定給定數(shù)據(jù)集的最佳組件數(shù)的自然方法诅愚。生成模型本質(zhì)上是數(shù)據(jù)集的概率分布寒锚，因此我們可以簡單地評估模型下數(shù)據(jù)的可能性，使用交叉驗證來避免過度擬合违孝。校正過度擬合的另一種方法是使用一些分析標(biāo)準(zhǔn)來調(diào)整模型可能性刹前，例如Akaike information criterion (AIC) 或 Bayesian information criterion (BIC)。Scikit-Learn的GMM估計器實際上包含計算這兩者的內(nèi)置方法雌桑，因此在這種方法上操作非常容易喇喉。
讓我們看看在moon數(shù)據(jù)集中，使用AIC和BIC函數(shù)確定GMM組件數(shù)量：

n_components = np.arange(1, 21)
models = [GMM(n, covariance_type='full', random_state=0).fit(Xmoon)
          for n in n_components]

plt.plot(n_components, [m.bic(Xmoon) for m in models], label='BIC')
plt.plot(n_components, [m.aic(Xmoon) for m in models], label='AIC')
plt.legend(loc='best')
plt.xlabel('n_components');

最佳的聚類數(shù)目是使得AIC或BIC最小化的值校坑，具體取決于我們希望使用的近似值拣技。 AIC告訴我們，我們上面選擇的16個組件可能太多了：大約8-12個組件可能是更好的選擇耍目。與此類問題一樣膏斤，BIC建議使用更簡單的模型。
注意重點：這個組件數(shù)量的選擇衡量GMM作為密度估算器的效果邪驮，而不是它作為聚類算法的效果莫辨。我鼓勵您將GMM主要視為密度估算器，并且只有在簡單數(shù)據(jù)集中保證時才將其用于聚類毅访。

例子：GMM生成新數(shù)據(jù)

我們剛剛看到了一個使用GMM作為數(shù)據(jù)生成模型的簡單示例沮榜，以便根據(jù)輸入數(shù)據(jù)定義的分布創(chuàng)建新樣本。在這里喻粹，我們將運行這個想法蟆融，并從我們以前使用過的標(biāo)準(zhǔn)數(shù)字語料庫中生成新的手寫數(shù)字。
首先磷斧，讓我們使用Scikit-Learn的數(shù)據(jù)工具加載數(shù)字?jǐn)?shù)據(jù)：

from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape

(1797, 64)

接下來讓我們繪制前100個振愿，以準(zhǔn)確回憶我們正在看的內(nèi)容：

def plot_digits(data):
    fig, ax = plt.subplots(10, 10, figsize=(8, 8),
                           subplot_kw=dict(xticks=[], yticks=[]))
    fig.subplots_adjust(hspace=0.05, wspace=0.05)
    for i, axi in enumerate(ax.flat):
        im = axi.imshow(data[i].reshape(8, 8), cmap='binary')
        im.set_clim(0, 16)
plot_digits(digits.data)

我們有64個維度的近1,800位數(shù)字捷犹，我們可以在這些位置上構(gòu)建GMM以產(chǎn)生更多冕末。 GMM可能難以在如此高維空間中收斂，因此我們將從數(shù)據(jù)上的可逆維數(shù)減少算法開始档桃。在這里，我們將使用一個簡單的PCA蔑舞，要求它保留99％的預(yù)測數(shù)據(jù)方差：

from sklearn.decomposition import PCA
pca = PCA(0.99, whiten=True)
data = pca.fit_transform(digits.data)
data.shape

(1797, 41)

結(jié)果是41個維度，減少了近1/3攻询，幾乎沒有信息丟失。根據(jù)這些預(yù)測數(shù)據(jù)钧栖，讓我們使用AIC來計算我們應(yīng)該使用的GMM組件的數(shù)量：

n_components = np.arange(50, 210, 10)
models = [GMM(n, covariance_type='full', random_state=0)
          for n in n_components]
aics = [model.fit(data).aic(data) for model in models]
plt.plot(n_components, aics);

似乎大約110個components最小化了AIC婆翔；我們將使用這個模型。我們迅速將其與數(shù)據(jù)擬合并確保它已收斂合：

gmm = GMM(110, covariance_type='full', random_state=0)
gmm.fit(data)
print(gmm.converged_)

True

現(xiàn)在我們可以使用GMM作為生成模型在這個41維投影空間內(nèi)繪制100個新點的樣本：

data_new = gmm.sample(100, random_state=0)
data_new.shape

(100, 41)

最后啃奴，我們可以使用PCA對象的逆變換來構(gòu)造新的數(shù)字：

digits_new = pca.inverse_transform(data_new)
plot_digits(digits_new)

大部分結(jié)果看起來像數(shù)據(jù)集中合理的數(shù)字潭陪！
考慮一下我們在這里做了什么：給定一個手寫數(shù)字的樣本，我們已經(jīng)模擬了數(shù)據(jù)的分布最蕾，這樣我們就可以從數(shù)據(jù)中生成全新的數(shù)字樣本：這些是“手寫數(shù)字”依溯，不是單獨的出現(xiàn)在原始數(shù)據(jù)集中，而是捕獲混合模型建模的輸入數(shù)據(jù)的一般特征揖膜。這種數(shù)字生成模型可以證明作為貝葉斯生成分類器的一個組成部分非常有用誓沸，我們將在下一節(jié)中看到。