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

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

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

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

GMM的動(dòng)機(jī)：k-means的弱點(diǎn)

我們看下k-means的缺陷，思考下如何提高聚類(lèi)模型。正如上一節(jié)所示，給定簡(jiǎn)單，易于分類(lèi)的數(shù)據(jù)，k-means能找到合適的聚類(lèi)結(jié)果。
舉例而言，假設(shè)我們有些簡(jiǎn)單的數(shù)據(jù)點(diǎn)，k-means算法能以某種方式很快地將它們聚類(lèi)，跟我們?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');

從直觀的角度來(lái)看，我可能期望聚類(lèi)分配時(shí)，某些點(diǎn)比其他的更確定：舉例而言，中間兩個(gè)聚類(lèi)之間似乎存在非常輕微的重疊，這樣我們可能對(duì)這些數(shù)據(jù)點(diǎn)的分配沒(méi)有完全的信心。不幸的是，k-means模型沒(méi)有聚類(lèi)分配的概率或不確定性的內(nèi)在度量（盡管可能使用bootstrap 的方式來(lái)估計(jì)這種不確定性）。為此，我們必須考慮泛化這種模型。
k-means模型的一種理解思路是，它在每個(gè)類(lèi)蔟的中心放置了一個(gè)圈（或者，更高維度超球面)，其半徑由聚類(lèi)中最遠(yuǎn)的點(diǎn)確定。該半徑充當(dāng)訓(xùn)練集中聚類(lèi)分配的一個(gè)硬截?cái)啵喝魏稳ν獾臄?shù)據(jù)點(diǎn)不被視為該類(lèi)的成員。我們可以使用以下函數(shù)可視化這個(gè)聚類(lèi)模型：

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的一個(gè)重要發(fā)現(xiàn)，這些聚類(lèi)模式必須是圓形的。k-means沒(méi)有內(nèi)置的方法來(lái)計(jì)算橢圓形或橢圓形的簇。因此，舉例而言，假設(shè)我們將相同的數(shù)據(jù)點(diǎn)作變換，這種聚類(lèi)分配方式最終變得混亂：

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）試圖找到一個(gè)多維高斯概率分布的混合，以模擬任何輸入數(shù)據(jù)集。在最簡(jiǎn)單的情況下，GMM可用于以與k-means相同的方式聚類(lèi)。

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');

但因?yàn)镚MM包含概率模型，因此可以找到聚類(lèi)分配的概率方式 - 在Scikit-Learn中，通過(guò)調(diào)用predict_proba方法實(shí)現(xiàn)。它將返回一個(gè)大小為[n_samples, n_clusters]的矩陣，用于衡量每個(gè)點(diǎn)屬于給定類(lèi)別的概率：

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.   ]]

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

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

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

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

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)

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

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

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

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

這清楚地表明GMM解決了以前遇到的k-means的兩個(gè)主要實(shí)際問(wèn)題。

選擇協(xié)方差類(lèi)型

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

GMM 作為密度估計(jì)

盡管GMM通常被歸類(lèi)為聚類(lèi)算法，但從根本上說(shuō)它是一種密度估算算法。也就是說(shuō)，GMM適合某些數(shù)據(jù)的結(jié)果在技術(shù)上不是聚類(lèi)模型，而是描述數(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]);

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

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

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

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

這里，16個(gè)高斯分布的混合不是為了找到分離的數(shù)據(jù)簇，而是為了對(duì)輸入數(shù)據(jù)的整體分布進(jìn)行建模。這是分布的一個(gè)生成模型，這意味著GMM為我們提供了生成與我們的輸入類(lèi)似分布的新隨機(jī)數(shù)據(jù)的方法。例如，以下是從這個(gè)16分量GMM擬合到我們?cè)紨?shù)據(jù)的400個(gè)新點(diǎn)：

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

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

多少components？

GMM是一種生成模型這一事實(shí)為我們提供了一種確定給定數(shù)據(jù)集的最佳組件數(shù)的自然方法。生成模型本質(zhì)上是數(shù)據(jù)集的概率分布，因此我們可以簡(jiǎn)單地評(píng)估模型下數(shù)據(jù)的可能性，使用交叉驗(yàn)證來(lái)避免過(guò)度擬合。校正過(guò)度擬合的另一種方法是使用一些分析標(biāo)準(zhǔn)來(lái)調(diào)整模型可能性，例如Akaike information criterion (AIC) 或 Bayesian information criterion (BIC)。Scikit-Learn的GMM估計(jì)器實(shí)際上包含計(jì)算這兩者的內(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');

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

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

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

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

(1797, 64)

接下來(lái)讓我們繪制前100個(gè)，以準(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個(gè)維度的近1,800位數(shù)字，我們可以在這些位置上構(gòu)建GMM以產(chǎn)生更多。 GMM可能難以在如此高維空間中收斂，因此我們將從數(shù)據(jù)上的可逆維數(shù)減少算法開(kāi)始。在這里，我們將使用一個(gè)簡(jiǎn)單的PCA，要求它保留99％的預(yù)測(cè)數(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個(gè)維度，減少了近1/3，幾乎沒(méi)有信息丟失。根據(jù)這些預(yù)測(cè)數(shù)據(jù)，讓我們使用AIC來(lái)計(jì)算我們應(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個(gè)components最小化了AIC；我們將使用這個(gè)模型。我們迅速將其與數(shù)據(jù)擬合并確保它已收斂合：

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

True

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

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

(100, 41)

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

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

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