高斯混合模型聚类_大咖聊技术|高斯混合模型

MathWorks

或许你知道如今企业要想在商业上拔得头筹，已经离不开AI的协助？

或许你听说过可以自动从数据中得到知识辅助决策的强大机器学习算法？

机器学习可以分成三个子领域：监督学习，无监督学习和强化学习。无监督学习可以看成是“没有老师情况下的学习”，因为只有数据本身，没有监督学习里的标签，也没有强化学习里的反馈。

这里我们介绍一种聚类方法，高斯混合模型(Gaussian mixture model)。(这里要和mixed model区分开来，mixed model是一种统计模型，主要用来处理重复测量，或者有群组效应的数据，可以认为是一种监督学习方法)。高斯混合模型是一种基于模型的聚类方法，它是混合模型的一种特例，因为它只考虑高斯分布的混合。它和另一种常见的聚类方法K-Means有很多相似性。从某种程度上来说，K-Means是高斯混合模型的一种特殊情况。能使用K-Means的问题就能使用高斯混合模型，比如图像压缩中使用的矢量量化(vector quantization，参见Elements of statistical learning, Chap 14)。除此之外，高斯混合模型还可以用来做异常检测。

模 型 定 义

具体模型如下：我们有一批数据{x_i},i=1,…,n,其中每个数据都是独立服从一个混合高斯分布。即

其中的概率密度函数可写成

这里的ϕ_k是一个高斯分布，有自己的参数μ_k和σ_k²。同时，x_i来自每个高斯分布的概率为π_k，满足

模型的参数有

满足

同时还有高斯分布的参数

此外还有一个超参数K。这里我们给出的记号是一维数据，但是推广到高维空间也并不难，只需要

其中⃗μ_k是p维空间向量，Σ_k是p×p半正定矩阵。在后续的论述中，我们会先使用一维数据为例子，然后再推广到高维空间。

因此，对于数据，我们可以写出对数似然函数

算 法 细 节

对于高斯混合模型来说，我们要估计的参数有：

·π_k,k∈{1,…,K}：各高斯组分比例π_k>0,∑_kπ_k=1

·μ_k,k∈{1,…,K}：各高斯组分均值

·σ_k²,k∈{1,…,K}：各高斯组分方差

我们求解方法是最大似然估计(maximum (log-)likelihood estimation, MLE)。对数似然函数是

对于高斯混合模型求解，比较好的方法是Expectation Maximization(EM)算法。对于每个观察值x_i来说，我们并不知道它属于哪个组分，为此，我们引入一个隐变量，组分指示器(indicator)z_ik。它的取值如下

当我们能够观测到它的值的时候，我们有全数据似然函数

相应的对数似然函数是

这时候参数估计非常简单。对于组分的均值和方差来说，就是简单的单高斯分布参数估计。对于组分比例来说，就是

但是实际情况是我们不知道z_ik的值，因此我们必须使用EM算法。EM算法大致步骤如下：

1.    初始化参数

2.    基于当前参数估计θ^(t)迭代以下两步直至收敛：

a.    E-step：计算在当前参数下隐变量z_ik的期望(expectation)，

b.    M-step：最大化对数似然函数在给定隐变量期望下的参数，并更新。

对于高斯混合模型来说，下面给出参数更新具体计算步骤。计算隐变量期望的表达式如下：

对于更新参数，我们有

对其中的参数求偏导并令其为零，并考虑约束条件∑_kπ_k=1可解得

上面的式子要推广到高维空间也很直接，只要保证相应的矩阵和矢量形状正确，不难写出正确的更新步骤。

算 法 实 现

下面我们用一维数据来实现高斯混合模型的EM算法，真实的参数是π₁=π₂=0.5，μ₁=−2，μ₂=2，σ₁²=σ₂²=1。

set.seed(1234)

x 2000, -2), rnorm(2000, 2))

hist(x, breaks=500,freq = FALSE)

对于初始化，我们采用一个十分简单的算法

pi1 0.5

pi2 0.5

mu1 -0.5

mu2 0.5

sd1

sd2

下面我们定义E-step。因为高斯分布很容易得到接近0的概率分布，为了解决这个数值稳定性问题，我们采用log-sum trick。

mat_x

  Estep function(param_list) {

    pi1

    pi2

    mu1

    mu2

    var1

    var2

    # compute aik

    mu_vec

    var_vec

    pi_vec

    aik

    aik 2/(2*var_vec)

    aik 0.5*log(var_vec)

    aik

    aik

    aik 1, max)

    gammaik

    gammaik 1, sum)

  }

我们现在可以定义M-step

Mstep function(gammaik) {

  nc

  pis

  pi1 1]

  pi2 2]

  mus

  mu1 1]

  mu2 2]

  vars 2 * gammaik) / nc

  var1 1]

  var2 2]

return(list(pi1=pi1, pi2=pi2, mu1=mu1, mu2=mu2, var1=var1,var2=var2))

}

有了这两个函数，我们开始估计模拟数据的参数了

N

param_list

  pi1=pi1, pi2=pi2,

  mu1=mu1, mu2=mu2,

  var1=sd1**2, var2=sd2**2)

for (i in1:100) {

  gammaik

  param_list

}

param_list

## $pi1

## [1] 0.499523

##

## $pi2

## [1] 0.500477

##

## $mu1

## [1] -2.0045

##

## $mu2

## [1] 2.004301

##

## $var1

## [1] 0.987977

##

## $var2

## [1] 1.023447

从这里我们可以看出，对于参数的估计比较准确。高斯混合模型还可以给出每个数据属于各个组分的概率Pr(z_ik=1|x,θ)，我们可以用这个来对数据进行聚类。

小 结

在本篇文章里，我们利用一维数据介绍了高斯混合模型的一种常用求解方法：EM算法，用代码展示了具体的计算步骤。

由于本文是介绍算法目的，代码更注重的是可读性，对于质量并没有过多追求。因此有很多可以改进的地方，比如说在E-step和M-step两个函数中都使用了全局变量，也没有收敛性的检查。在实际生产代码中，这些都是需要改进的地方。

读完这篇分享，数据科学爱好者们，是不是感觉到机器学习是一个很有意思的领域，想要进一步了解呢？

- END -

本期大咖员工

Dr. Liu

PhD, Physics, Université Paris Sud

Master, Statistics, Université Paris Sud

Master, Physics, Université Paris Sud

Bachelor, Physics, Peking University