高斯混合模型（前面本文中将应用他的缩写 GMM）听起来很简单，其实他的工作原理和 KMeans 十分类似，你甚至能够认为它是 KMeans 的概率版本。 这种概率特色使 GMM 能够利用于 KMeans 无奈解决的许多简单问题。

因为KMeans的限度很多，比方： 它假如簇是球形的并且大小雷同，这在大多数事实世界的场景中是有效的。并且它是硬聚类办法，这意味着每个数据点都调配给一个集群，这也是不事实的。

在本文中，咱们将依据下面的内容来介绍 KMeans 的一个代替计划之一，高斯混合模型。

从概念上解释：高斯混合模型就是用高斯概率密度函数（正态分布曲线）准确地量化事物，它是一个将事物合成为若干的基于高斯概率密度函数（正态分布曲线）造成的模型。

高斯混合模型 (GMM) 算法的工作原理

正如后面提到的，能够将 GMM 称为 概率的KMeans，这是因为 KMeans 和 GMM 的终点和训练过程是雷同的。 然而，KMeans 应用基于间隔的办法，而 GMM 应用概率办法。 GMM 中有一个次要假如：数据集由多个高斯分布组成，换句话说，GMM 模型能够看作是由 K 个单高斯模型组合而成的模型,这 K 个子模型是混合模型的隐变量(Hidden variable)

上述散布通常称为多模型散布。 每个峰代表咱们数据集中不同的高斯分布或聚类。 咱们肉眼能够看到这些散布，然而应用公式如何预计这些散布呢？

在解释这个问题之前，咱们先创立一些高斯分布。这里咱们生成的是多元正态分布； 它是单变量正态分布的更高维扩大。

让咱们定义数据点的均值和协方差。 应用均值和协方差，咱们能够生成如下散布。

 # Set the mean and covariance mean1= [0, 0] mean2= [2, 0] cov1= [[1, .7], [.7, 1]] cov2= [[.5, .4], [.4, .5]]  # Generate data from the mean and covariance data1=np.random.multivariate_normal(mean1, cov1, size=1000) data2=np.random.multivariate_normal(mean2, cov2, size=1000)

可视化一下生成的数据

 plt.figure(figsize=(10,6))  plt.scatter(data1[:,0],data1[:,1]) plt.scatter(data2[:,0],data2[:,1])  sns.kdeplot(data1[:, 0], data1[:, 1], levels=20, linewidth=10, color='k', alpha=0.2) sns.kdeplot(data2[:, 0], data2[:, 1], levels=20, linewidth=10, color='k', alpha=0.2)  plt.grid(False) plt.show()

咱们下面所做的工作是：应用均值和协方差矩阵生成了随机高斯分布。 而 GMM 要做正好与这个相同，也就是找到一个散布的均值和协方差，那么怎么做呢？

工作过程大抵如下：

为给定的数据集确定聚类的数量（这里咱们能够应用畛域常识或其余办法，例如 BIC/AIC）。 依据咱们下面的参数，有 1000 个数据点，和两个簇2。

初始化每个簇的均值、协方差和权重参数。

应用冀望最大化算法执行以下操作：

  • 冀望步骤（E-step）：计算每个数据点属于每个散布的概率，而后应用参数的以后预计评估似然函数
  • 最大化步骤（M-step）：更新之前的均值、协方差和权重参数，这样最大化E步骤中找到的预期似然
  • 反复这些步骤，直到模型收敛。

以上是GMM 算法的非数学的通俗化的解释。

GMM 数学原理

有了下面的艰深解释，咱们开始进入正题，下面的解释能够看到GMM 的外围在于上一节中形容的冀望最大化 (EM) 算法。

在解释之前，咱们先演示一下 EM 算法在 GMM 中的利用。

1、初始化均值、协方差和权重参数

  • mean (): 随机初始化
  • 协方差 ()：随机初始化
  • 权重（混合系数）()：每个类的分数是指特定数据点属于每个类的可能性。 一开始，这对所有簇都是平等的。 假如咱们用三个重量拟合 GMM，那么每个组件的权重参数可能设置为 1/3，这样概率分布为 (1/3, 1/3, 1/3)。

2、冀望步骤（E-step）

对于每个数据点，应用以下等式计算数据点属于簇 () 的概率。 这里的k是我散布（簇）数。

下面的_是高斯分布c的混合系数（有时称为权重），它在上一阶段被初始化，（|，）形容了高斯分布的概率密度函数（PDF），均值为和 对于数据点 x 的协方差 ； 所以能够有如下示意。

E-step 应用模型参数的以后估计值计算概率。 这些概率通常称为高斯分布的“responsibilities”。 它们由变量 r_ic 示意，其中 i 是数据点的索引，c 是高斯分布的索引。 responsibilities掂量第 c 个高斯分布对生成第 i 个数据点“负责”的水平。 这里应用条件概率，更具体地说是贝叶斯定理。

一个简略的例子。 假如咱们有 100 个数据点，须要将它们聚类成两组。 咱们能够这样写 r_ic(i=20,c=1) 。 其中 i 代表数据点的索引，c 代表咱们正在思考的簇的索引。

不要遗记在开始时，_ 会初始化为等值。 在咱们的例子中，_1 = _2 = 1/2。

E-step 的后果是混合模型中每个数据点和每个高斯分布的一组responsibilities。 这些responsibilities会在 M-step更新模型参数的预计。

3、最大化步骤（M-step）

算法应用高斯分布的responsibilities（在 E-step中计算的）来更新模型参数的估计值。

M-step更新参数的预计如下：

下面的公式 4 更新 c（混合系数），公式 5 更新 c，公式6 更新 c。更新后的的预计会在下一个 E-step 中用于计算数据点的新responsibilities。

GMM 将将反复这个过程直到算法收敛，通常在模型参数从一次迭代到下一次迭代没有显着变动时就被认为收敛了。

让咱们将下面的过程整顿成一个简略的流程图，这样能够更好的了解：

数学原理完了，上面该开始应用 Python 从头开始实现 GMM了

应用 Python 从头开始实现 GMM

了解了数学原理，GMM的代码也不简单，基本上下面的每一个公式应用1-2行就能够实现

首先，创立一个试验数据集，咱们将为一维数据集实现 GMM，因为这个比较简单

 importnumpyasnp  n_samples=100 mu1, sigma1=-5, 1.2 mu2, sigma2=5, 1.8 mu3, sigma3=0, 1.6  x1=np.random.normal(loc=mu1, scale=np.sqrt(sigma1), size=n_samples) x2=np.random.normal(loc=mu2, scale=np.sqrt(sigma2), size=n_samples) x3=np.random.normal(loc=mu3, scale=np.sqrt(sigma3), size=n_samples)  X=np.concatenate((x1,x2,x3))

上面是可视化的函数封装

 fromscipy.statsimportnorm  defplot_pdf(mu,sigma,label,alpha=0.5,linestyle='k--',density=True):     """     Plot 1-D data and its PDF curve.      """     # Compute the mean and standard deviation of the data      # Plot the data          X=norm.rvs(mu, sigma, size=1000)          plt.hist(X, bins=50, density=density, alpha=alpha,label=label)      # Plot the PDF     x=np.linspace(X.min(), X.max(), 1000)     y=norm.pdf(x, mu, sigma)     plt.plot(x, y, linestyle)

并绘制生成的数据如下。 这里没有绘制数据自身，而是绘制了每个样本的概率密度。

 plot_pdf(mu1,sigma1,label=r"$\mu={} \ ; \ \sigma={}$".format(mu1,sigma1)) plot_pdf(mu2,sigma2,label=r"$\mu={} \ ; \ \sigma={}$".format(mu2,sigma2)) plot_pdf(mu3,sigma3,label=r"$\mu={} \ ; \ \sigma={}$".format(mu3,sigma3)) plt.legend() plt.show()

上面开始依据下面的公式实现代码：

1、初始化

 defrandom_init(n_compenents):          """Initialize means, weights and variance randomly        and plot the initialization     """          pi=np.ones((n_compenents)) /n_compenents     means=np.random.choice(X, n_compenents)     variances=np.random.random_sample(size=n_compenents)      plot_pdf(means[0],variances[0],'Random Init 01')     plot_pdf(means[1],variances[1],'Random Init 02')     plot_pdf(means[2],variances[2],'Random Init 03')          plt.legend()     plt.show()          returnmeans,variances,pi

2、E step

 defstep_expectation(X,n_components,means,variances):     """E Step          Parameters     ----------     X : array-like, shape (n_samples,)         The data.     n_components : int         The number of clusters     means : array-like, shape (n_components,)         The means of each mixture component.     variances : array-like, shape (n_components,)         The variances of each mixture component.              Returns     -------     weights : array-like, shape (n_components,n_samples)     """     weights=np.zeros((n_components,len(X)))     forjinrange(n_components):         weights[j,:] =norm(loc=means[j],scale=np.sqrt(variances[j])).pdf(X)     returnweights

3、M step

 defstep_maximization(X,weights,means,variances,n_compenents,pi):     """M Step          Parameters     ----------     X : array-like, shape (n_samples,)         The data.     weights : array-like, shape (n_components,n_samples)         initilized weights array     means : array-like, shape (n_components,)         The means of each mixture component.     variances : array-like, shape (n_components,)         The variances of each mixture component.     n_components : int         The number of clusters     pi: array-like (n_components,)         mixture component weights              Returns     -------     means : array-like, shape (n_components,)         The means of each mixture component.     variances : array-like, shape (n_components,)         The variances of each mixture component.     """     r= []     forjinrange(n_compenents):            r.append((weights[j] *pi[j]) / (np.sum([weights[i] *pi[i] foriinrange(n_compenents)], axis=0)))          #5th equation above         means[j] =np.sum(r[j] *X) / (np.sum(r[j]))                  #6th equation above         variances[j] =np.sum(r[j] *np.square(X-means[j])) / (np.sum(r[j]))                  #4th equation above         pi[j] =np.mean(r[j])      returnvariances,means,pi

而后就是训练的循环

 deftrain_gmm(data,n_compenents=3,n_steps=50, plot_intermediate_steps_flag=True):     """ Training step of the GMM model          Parameters     ----------     data : array-like, shape (n_samples,)         The data.     n_components : int         The number of clusters     n_steps: int         number of iterations to run     """          #intilize model parameters at the start     means,variances,pi=random_init(n_compenents)      forstepinrange(n_steps):         #perform E step         weights=step_expectation(data,n_compenents,means,variances)         #perform M step         variances,means,pi=step_maximization(X, weights, means, variances, n_compenents, pi)      plot_pdf(means,variances)

这样就实现了，让咱们看看 GMM 体现如何。

在上图中红色虚线示意原始散布，而其余色彩示意学习散布。 第 30 次迭代后，能够看到的模型在这个生成的数据集上体现良好。

这里只是为了解释GMM的概念进行的Python实现，在理论用例中请不要间接应用，请应用scikit-learn提供的GMM，因为它比咱们这个手写的要快多了，具体的对象名是

sklearn.mixture.GaussianMixture

,它的罕用参数如下：

  • tol：定义模型的进行规范。 当下限均匀增益低于 tol 参数时，EM 迭代将进行。
  • init_params：用于初始化权重的办法

总结

本文对高斯混合模型进行全面的介绍，心愿浏览完本文后你对 GMM 可能有一个具体的理解，GMM 的一个常见问题是它不能很好地扩大到大型数据集。如果你想本人尝试，本文的残缺代码在这里：https://avoid.overfit.cn/post/c9ca4de6d42d4d67b80cc789e921c636

作者：Ransaka Ravihara