python线性模型评价（高斯混合模型的解释及Python实现）

在机器学习领域，我们可以区分为两个主要领域:监督学习和无监督学习。两者的主要区别在于数据的性质以及处理数据的方法。聚类是一种无监督学习问题，我们希望在数据集中找到具有某些共同特征的点的聚类。假设我们有这样一个机器学习数据集:

我们的工作是找到一组看起来很接近的点。在这种情况下，我们可以清楚地识别出两组点，我们将分别把它们标记成蓝色和红色:

请注意，我们现在引入了一些额外的符号。这里，μ1和μ2是每个聚类的质心。一种流行的聚类算法为K-means，它将遵循迭代方法来更新每个聚类的参数。更具体地说，它将做的是计算每个聚类的均值（或质心），然后计算它们与每个数据点的距离。然后将后者标记为由其最接近的质心识别的聚类的一部分。重复此过程，直到满足某些收敛标准，例如，当我们看不到聚类分配中的进一步更改时。

K-means的一个重要特点是它是一种硬聚类方法，这意味着它将每个点关联到一个且只有一个聚类。这种方法的一个局限性在于，当有一些点过于靠近给定聚类的边界时，这些点可能会被错误地标记。那么，如何使用软聚类而不是硬聚类呢?这正是高斯混合模型（或简称为GMM）试图做的。现在让我们进一步讨论这个方法。

定义

高斯混合函数是由几个高斯函数组成的函数，每个由k∈{1，…，k}标识，其中k是我们数据集的聚类数。每个高斯函数由以下参数组成:

• 均值μₖ定义它的中心。
• 协方差Σₖ定义其宽度。这相当于多变量情况下椭球的尺寸。
• 混合概率πₖ定义高斯函数的大小。

现在让我们以图形方式说明这些参数：

在这里，我们可以看到有三个高斯函数，因此K = 3。每个高斯解释了三个可用聚类中包含的数据。每个高斯由平均μₖ，协方差Σₖ和混合系数πₖ定义。混合系数本身是概率，必须满足以下条件:

现在我们如何确定这些参数的最佳值呢？为了实现这一点，我们必须确保每个高斯函数都适合属于每个聚类的数据点。这正是极大似然的作用。

通常，高斯密度函数由下式给出：

其中x代表我们的数据点，D是每个数据点的维数。μ和Σ分别是均值和协方差。如果我们有一个由N = 1000个三维点（D = 3）组成的数据集，那么x将是一个1000×3矩阵。μ将是1×3向量，Σ将是3×3矩阵。我们还会发现得到这个等式的对数很有用，该等式由下式给出：

如果我们推导出这个关于均值和协方差的方程，然后把它等于零，那么我们就能找到这些参数的最优值，这些解对应于这个设置的最大似然估计(MLE)。但是，因为我们要处理的不只是一个高斯函数，而是很多高斯函数，当我们找到整个混合函数的参数时，事情就会变得有点复杂。

初步推导

首先，让我们假设一个数据点的概率X ₙ来自高斯ķ。我们可以表达为：

“ 给一个数据点X ₙ，它来自高斯k的概率是多少?”在这种情况下，zₖ是一个潜在变量，只有两个可能的值。X ₙ高斯来到ķ，否则为0。我们实际上并没有看到这个z变量，但是知道它的发生概率将有助于我们确定高斯混合参数，我们将在后面讨论。

我们可以这样说:

这意味着从高斯k中观察到一个点的总概率实际上等于这个高斯函数的混合系数。这是有道理的，因为高斯函数越大，我们期望的概率就越高。现在z为所有可能的zₖ的集合，因此

我们事先知道每个zₖ独立于其他z，并且当k等于点来自的聚类时它们只能取值1 。因此：

现在，如果我们的数据来自高斯k，那么要找出观察数据的概率呢?事实证明它实际上是高斯函数本身！按照我们用来定义p(z)的相同逻辑，我们可以声明:

现在你可能会问，我们为什么要做这些呢？还记得我们最初的目的是，以确定哪些概率zₖ（给我们的观察X ₙ）？事实证明，我们刚刚导出的方程以及贝叶斯规则将帮助我们确定这个概率。根据概率的乘积法则，我们知道

因此

这是一个定义高斯混合的方程，你可以清楚地看到它依赖于我们之前提到的所有参数!为了确定这些值的最优值，我们需要确定模型的最大似然值。我们能找到的可能性作为所有观测xₙ的联合概率,定义如下：

就像我们对原始高斯密度函数所做的那样，让我们​​将对数应用到等式的每一边：

现在为了找到高斯混合的最佳参数，我们所要做的就是得到关于πₖ，μₖ和Σₖ的这个方程式。我们这里有一个问题。计算此表达式的导数然后求解参数将非常困难！

我们需要使用迭代方法来估计参数。但首先，请记住我们应该找到p（zₖ = 1 | X ₙ），因为我们已经准备好了所有的东西来定义这个概率。

从贝叶斯规则来看，我们知道这一点

从我们以前的推导中，我们知道:

所以让我们替换这些p (zₖ= 1 | xₙ):

这就是我们一直在寻找的东西！接下来，我们将继续讨论一种有助于我们轻松确定高斯混合参数的方法。

在Python中实现

让我们首先导入Python库:

import imageio import matplotlib.animation as ani import matplotlib.cm as cmx import matplotlib.colors as colors import matplotlib.pyplot as plt import numpy as np from matplotlib.patches import Ellipse from PIL import Image from sklearn import datasets from sklearn.cluster import KMeans

我们将使用Iris机器学习数据集，我们可以使用sklearn提供的load_iris函数轻松获取:

iris = datasets.load_iris() X = iris.data

现在我们来实现高斯密度函数。

def gaussian(X, mu, cov): n = X.shape[1] diff = (X - mu).T return np.diagonal(1 / ((2 * np.pi) ** (n / 2) * np.linalg.det(cov) ** 0.5) * np.exp(-0.5 * np.dot(np.dot(diff.T, np.linalg.inv(cov)), diff))).reshape(-1, 1) x0 = np.array([[0.05, 1.413, 0.212], [0.85, -0.3, 1.11], [11.1, 0.4, 1.5], [0.27, 0.12, 1.44], [88, 12.33, 1.44]]) mu = np.mean(x0, axis=0) cov = np.dot((x0 - mu).T, x0 - mu) / (x0.shape[0] - 1) y = gaussian(x0, mu=mu, cov=cov) y

Step 1

这是GMM的初始化步骤。此时，我们必须初始化参数πₖ、μₖ和Σₖ。在本例中，我们将使用KMeans的结果作为μₖ的初始值，将πₖ设置为聚类数量1，并将Σₖ设置为单位矩阵。我们也可以使用随机数，但是使用一个合理的初始化过程将帮助算法获得更好的结果。

def initialize_clusters(X, n_clusters): clusters = [] idx = np.arange(X.shape[0]) # We use the KMeans centroids to initialise the GMM kmeans = KMeans().fit(X) mu_k = kmeans.cluster_centers_ for i in range(n_clusters): clusters.append({ 'pi_k': 1.0 / n_clusters, 'mu_k': mu_k[i], 'cov_k': np.identity(X.shape[1], dtype=np.float64) }) return clusters

Step 2 (Expectation step)

我们可以通过以下表达式来实现:

为了方便起见，我们只是将分母计算为分子中所有项的和，然后将其分配给一个名为total的变量

def expectation_step(X, clusters): totals = np.zeros((X.shape[0], 1), dtype=np.float64) for cluster in clusters: pi_k = cluster['pi_k'] mu_k = cluster['mu_k'] cov_k = cluster['cov_k'] gamma_nk = (pi_k * gaussian(X, mu_k, cov_k)).astype(np.float64) for i in range(X.shape[0]): totals[i] = gamma_nk[i] cluster['gamma_nk'] = gamma_nk cluster['totals'] = totals for cluster in clusters: cluster['gamma_nk'] /= cluster['totals']

Step 3 (Maximization step):

现在我们来实现最大化步骤。我们可以简单地定义:

然后我们可以通过以下方法计算修正后的参数:

注意:为了计算协方差，我们定义了一个包含(xn-μk)T的辅助变量diff。

def maximization_step(X, clusters): N = float(X.shape[0]) for cluster in clusters: gamma_nk = cluster['gamma_nk'] cov_k = np.zeros((X.shape[1], X.shape[1])) N_k = np.sum(gamma_nk, axis=0) pi_k = N_k / N mu_k = np.sum(gamma_nk * X, axis=0) / N_k for j in range(X.shape[0]): diff = (X[j] - mu_k).reshape(-1, 1) cov_k = gamma_nk[j] * np.dot(diff, diff.T) cov_k /= N_k cluster['pi_k'] = pi_k cluster['mu_k'] = mu_k cluster['cov_k'] = cov_k

我们现在确定模型的对数似然。由下式给出

然而，第二个求和已经在expectation_step函数中计算过，并且可以在total变量中使用。

def get_likelihood(X, clusters): likelihood = [] sample_likelihoods = np.log(np.array([cluster['totals'] for cluster in clusters])) return np.sum(sample_likelihoods), sample_likelihoods

最后，让我们把所有的东西放在一起!首先，我们将使用initialise_clusters函数初始化参数，然后执行几个expectation-maximization steps。在本例中，我们将训练过程的迭代次数设置为固定的n_epochs次数。我这样做是为了以后生成对数似然图。

def train_gmm(X, n_clusters, n_epochs): clusters = initialize_clusters(X, n_clusters) likelihoods = np.zeros((n_epochs, )) scores = np.zeros((X.shape[0], n_clusters)) history = [] for i in range(n_epochs): clusters_snapshot = [] # This is just for our later use in the graphs for cluster in clusters: clusters_snapshot.append({ 'mu_k': cluster['mu_k'].copy(), 'cov_k': cluster['cov_k'].copy() }) history.append(clusters_snapshot) expectation_step(X, clusters) maximization_step(X, clusters) likelihood, sample_likelihoods = get_likelihood(X, clusters) likelihoods[i] = likelihood print('Epoch: ', i 1, 'Likelihood: ', likelihood) for i, cluster in enumerate(clusters): scores[:, i] = np.log(cluster['gamma_nk']).reshape(-1) return clusters, likelihoods, scores, sample_likelihoods, history

训练模型

n_clusters = 3 n_epochs = 50 clusters, likelihoods, scores, sample_likelihoods, history = train_gmm(X, n_clusters, n_epochs)

查看对数似然

plt.figure(figsize=(10, 10)) plt.title('Log-Likelihood') plt.plot(np.arange(1, n_epochs 1), likelihoods) plt.show()

我们知道创建一个图形来可视化我们的聚类以及高斯混合的参数。实际上，我们所做的是创建不同尺度的椭圆以便映射到每个高斯函数的坐标。

def create_cluster_animation(X, history, scores): fig, ax = plt.subplots(1, 1, figsize=(10, 10)) colorset = ['blue', 'red', 'black'] images = [] for j, clusters in enumerate(history): idx = 0 if j % 3 != 0: continue plt.cla() for cluster in clusters: mu = cluster['mu_k'] cov = cluster['cov_k'] eigenvalues, eigenvectors = np.linalg.eigh(cov) order = eigenvalues.argsort()[::-1] eigenvalues, eigenvectors = eigenvalues[order], eigenvectors[:, order] vx, vy = eigenvectors[:,0][0], eigenvectors[:,0][1] theta = np.arctan2(vy, vx) color = colors.to_rgba(colorset[idx]) for cov_factor in range(1, 4): ell = Ellipse(xy=mu, width=np.sqrt(eigenvalues[0]) * cov_factor * 2, height=np.sqrt(eigenvalues[1]) * cov_factor * 2, angle=np.degrees(theta), linewidth=2) ell.set_facecolor((color[0], color[1], color[2], 1.0 / (cov_factor * 4.5))) ax.add_artist(ell) ax.scatter(cluster['mu_k'][0], cluster['mu_k'][1], c=colorset[idx], s=1000, marker=' ') idx = 1 for i in range(X.shape[0]): ax.scatter(X[i, 0], X[i, 1], c=colorset[np.argmax(scores[i])], marker='o') fig.canvas.draw() image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8') image = image.reshape(fig.canvas.get_width_height()[::-1] (3,)) images.append(image) kwargs_write = {'fps':1.0, 'quantizer':'nq'} imageio.mimsave('./gmm.gif', images, fps=1) plt.show(Image.open('gmm.gif').convert('RGB')) create_cluster_animation(X, history, scores)

让我们知道我们的计算是否正确。在本例中，我们使用sklearn的GMM实现来检查参数和概率。

from sklearn.mixture import GaussianMixture gmm = GaussianMixture(n_components=n_clusters, max_iter=50).fit(X) gmm_scores = gmm.score_samples(X) print('Means by sklearn:\n', gmm.means_) print('Means by our implementation:\n', np.array([cluster['mu_k'].tolist() for cluster in clusters])) print('Scores by sklearn:\n', gmm_scores[0:20]) print('Scores by our implementation:\n', sample_likelihoods.reshape(-1)[0:20])

,