EM算法与混合高斯模型(GMM)

1. 预备知识

1.1 期望(Expectation)

在概率论和统计学中，一个离散型随机变量的期望值（或数学期望、或均值，亦简称期望，物理中称为期望值）是试验中每次可能结果的概率乘以其结果的总和。

采用形式化定义，设\(Y\)是随机变量\(X\)的函数，\(Y=g(X)\)(\(g\)是连续函数)，那么

1）\(X\)是离散型随机变量，它的分布律为\(P(X=x_k)=p_k, k=1,2,\cdots\)，若\(\sum_{k=1}^{\infty}g(x_k)p_k\)绝对收敛，则期望值计算为

\[E[Y] = E[g(X)] = \sum_{k=1}^{\infty}g(x_k)p_k \tag{1.1} \]

2）\(X\)是连续型随机变量，存在相应的概率密度函数\(f(x)\)，若积分\(\int_{-\infty}^{+\infty}g(x)f(x)dx\)绝对收敛，则期望值计算为

\[E[Y] = E[g(x)] = \int_{-\infty}^{+\infty}g(x)f(x)dx \tag{1.2} \]

1.2 Jensen不等式

首先，介绍凸函数(convex function)的概念，如果一个函数具有以下的性质：每条弦都位于函数图像或其上⽅（如图所⽰），那么我们说这个函数是凸函数。位于\(x = a\)到\(x = b\)之间的任何⼀个\(x\)值都可以写成\(\lambda{a}+(1-\lambda)b\)的形式，其中\(0 \le \lambda \le 1\)。弦上的对应点可以写成\(\lambda {f(a)} + (1 -\lambda)f(b)\)，函数的对应值为\(f(\lambda a + (1 -\lambda)b)\)。这样，凸函数的性质就可以表⽰为:

\[f(\lambda a + (1 -\lambda)b) \le \lambda {f(a)} + (1 -\lambda)f(b) \tag{1.3} \]

这等价于要求函数的⼆阶导数处处为正。凸函数的例⼦有\(xln x(x > 0)\)和\(x^2\) 。如果等号只在\(\lambda = 0\)和\(\lambda = 1\)处取得，我们就说这个函数是严格凸函数（strictly convex function）。如果⼀个函数具有相反的性质，即每条弦都位于函数图像或其下⽅，那么这个函数被称为凹函数（concave function）。对应地，也有严格凹函数（strictly concave function）的定义。如果\(f(x)\)是凸函数，那么\(-f(x)\)就是凹函数。

图1.1 凸函数f(x)的每条弦（蓝⾊表⽰）位于函数上或函数上⽅，函数⽤红⾊曲线表⽰

使⽤归纳法，我们可以根据公式\((1.3)\)证明凸函数\(f(x)\)满⾜

\[f(\sum_{i=1}^M \lambda_i x_i) \le \sum_{i=1}^M \lambda_i f(x_i) \tag{1.4} \]

其中，对于任意点集\(\{x_i\}\)，都有\(\lambda_i \ge 0\)且\(\sum_i \lambda_i=1\)。上述公式的结果被称为Jensen不等式(Jensen's inequality)。如果我们把\(\lambda_i\)看成取值为\(\{x_i\}\)的离散变量\(x\)的概率分布，那么上面的公式可以写成

\[f(E[x]) \le E[f(x)] \tag{1.5} \]

其中，\(E[\cdot]\)表示期望。对于连续变量，Jensen不等式的形式为

\[f(\int xp(x)dx) \le \int f(x)p(x)dx \tag{1.6} \]

如果\(f(x)\)是凹函数，则

\[f(E[x]) \ge E[f(x)] \tag{1.7} \]

在EM算法中，\(f(x)\)为\(log\)函数（凹函数），将用到此不等式，特此强调。

特别的，如果\(f(x)\)是严格凸函数，那么\(f(E[x]) = E[f(x)]\)当且仅当\(P(X=E[x])=1\)，也就是说\(X\)是一个常量。

1.3 极大似然估计

问题描述

假如我们需要调查学校的男生和女生的身高分布，我们抽取100个男生和100个女生，将他们按照性别划分为两组。然后，统计抽样得到100个男生的身高数据和100个女生的身高数据。如果我们知道他们的身高服从正态分布，但是这个分布的均值 \(\mu\)和方差\(\sigma^2\)是不知道的，这两个参数就是我们需要估计的。

问题：我们已知的条件有两个，1）样本服从的概率分布模型，2）随机抽取的一些样本。那么，我们怎么求解该模型的参数？

我们可以根据已知条件，通过极大似然估计，求出未知参数。极大似然估计(maximum likelihood estimation, MLE)是一种模型参数估计的常用方法。这一方法的使用情境常常是这样的，对于给定的模型以及已经观察到的样本值\(X=\{x_1,x_2,\cdots,x_n\}\)，依据已经得到的样本的观察值，来估计该模型的参数\(\theta\)。而其中蕴含的一个直观的想法是：在给定模型下，现在已经得到样本值\(\{x_1,x_2,\cdots,x_n\}\)，这表示取得这一观察值的概率比较大，而我们所估计出来的参数，正是为了使现有观察情况出现可能性最大。要特别说明的是，极大似然估计有一个很重要的假设，就是所使用的样本之间是满足独立同分布的。

用数学知识解决现实问题

根据极大似然估计方法的描述，我们将前面的问题数学化表示：样本集\(X=\{x_1,x_2,\cdots,x_n\}, n=100\)。概率密度是：\(p(x_i|\theta)\)抽到第\(i\)个男生身高的概率。由于100个样本之间独立同分布，所以同时抽到这100个男生的概率是它们各自概率的乘积，就是从分布是\(p(X|\theta)\)的总体样本中抽到这100个样本的概率，也就是样本集\(X\)中各个样本的联合概率，用下式表示：

\[L(\theta) = L(x_1,x_2,\cdots,x_n;\theta) = \prod_{i=1}^np(x_i;\theta), \theta \in \Theta \tag{1.8} \]

这个概率反映了在概率密度函数的参数是\(\theta\)时，得到\(X\)这组样本的概率。我们需要找到一个参数\(\theta\)，使得抽到\(X\)这组样本的概率最大，也就是说需要其对应的似然函数\(L(\theta)\)最大。满足条件的\(\theta\)叫做\(\theta\)的极大似然估计值，记为：

\[\hat{\theta} = argmaxL(\theta) \tag{1.9} \]

极大似然估计值的求解步骤

1）写出似然函数

\[L(\theta) = L(x_1,x_2,\cdots,x_n;\theta) = \prod_{i=1}^np(x_i;\theta), \theta \in \Theta \]

2）对似然函数取对数

\[l(\theta) = ln \, L(\theta) = ln \, \prod_{i=1}^n p(x_i;\theta) = \sum_{i=1}^nln \, p(x_i;\theta) \]

3）对上式求导，令导数为0，得到似然方程。

4）解似然方程，得到的参数值即为所求。

多数情况下，我们是根据已知条件来推算结果，而极大似然估计是已知结果，寻求使该结果出现的可能性最大的条件，以此作为估计值。

2. EM算法引入

假设有1枚硬币A，这枚硬币正面出现的概率是\(p\)。进行掷硬币试验：投掷硬币A，出现正面记作1，出现反面记作0；独立地重复\(n\)次。这是一个参数为\(p\)的伯努利分布：

\[p(x) = p^x(1-p)^{1-x} \]

如果给定一个结果序列\(n=10\)，如下

\[1 \quad 1\quad 0\quad 1\quad 1\quad 1\quad 0\quad 0\quad 1\quad 1 \]

\(p\)可以通过极大似然估计得到

\[p=\cfrac{\sum_{i=1}^nI(x_i=1)}{n} = \cfrac{7}{10} \]

其中，\(I(\cdot)\)是指示函数。

现在考虑更复杂的情况，假设有3枚硬币，分别记作A，B，C。这些硬币正面出现的概率分别是\(\pi\)，\(p\)和\(q\)。进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或硬币C，如果是正面选择硬币B，如果是反面选择硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复\(n\)次试验（这里，\(n=10\)），观察结果如下：

\[1 \quad 1\quad 0\quad 0\quad 1\quad 0\quad 1\quad 0\quad 1\quad 1 \]

这种情况下，如果我们知道了如下带有颜色的结果，红色表示由硬币B掷出，绿色表示由硬币C掷出，即我们知道每次硬币A掷出的结果：

\[\color{red}{1} \quad \color{red}{1}\quad \color{green}{0}\quad \color{green}{0}\quad \color{red}{1}\quad \color{green}{0}\quad \color{red}{1}\quad \color{red}{0}\quad \color{green}{1}\quad \color{red}{1} \]

我们可以很容易的估计出参数\(\pi\)，\(p\)和\(q\):

\[\pi = \cfrac{红色的个数}{总个数} \quad p=\cfrac{红色中1的个数}{红色的个数} \quad q=\cfrac{绿色中1的个数}{绿色的个数} \]

但是现在的问题是，我们并不知道每次硬币A掷出的结果，也就是我们并不能够看到如上的带有颜色的结果，这时怎么估计参数\(\pi\)，\(p\)和\(q\)呢？

对于这个问题，按照极大似然估计的思路，令\(\theta=(\pi,p,q)\)，我们希望求得参数\(\theta\)，使得

\[L(\theta) = arg\smash{\max_\theta}p(x;\theta) \]

但是因为此时我们无法观测到掷硬币A的结果，所以需要引入隐变量\(z\)，\(z \in \{0,1\}\)，若\(z=1\)表示使用硬币B，若\(z=0\)表示使用硬币C，因此

\[\begin{aligned} p(x;\theta) &= \sum_{z}p(x,z;\theta) \\ &= \sum_{z}p(x|z;\theta)p(z;\theta) \\ &= p(x|z=1;\theta)p(z=1;\theta) + p(x|z=0;\theta)p(z=0;\theta) \end{aligned}\]

其中，我们知道

\[p(x|z=1;\theta) = p^x(1-p)^{1-x} \]

\[p(x|z=0;\theta) = q^x(1-q)^{1-x} \]

\[p(z=1;\theta)=\pi \quad p(z=0;\theta)=1-\pi \]

所以有，

\[p(x;\theta) = {\color{red}{\pi p^x(1-p)^{1-x}}} + \color{green}{(1-\pi)q^x(1-q)^{1-x}} \]

我们用\(\mu\)表示在给定当前参数\(\theta=(\pi,p,q)\)情况下，每个试验结果\(x\)由硬币B产生的概率，则

\[\mu=\cfrac{\color{red}{\pi p^x(1-p)^{1-x}}}{{\color{red}{\pi p^x(1-p)^{1-x}}} + \color{green}{(1-\pi)q^x(1-q)^{1-x}}}, x\in \{0,1\} \]

给定一个观测到的结果序列\(x=\{x_1,x_2,\cdots，x_N\}\)，则有

\[\mu_i=\cfrac{\color{red}{\pi p^{x_i}(1-p)^{1-{x_i}}}}{{\color{red}{\pi p^{x_i}(1-p)^{1-{x_i}}}} + \color{green}{(1-\pi)q^{x_i}(1-q)^{1-{x_i}}}} \]

至此，需要发挥一点想象力，假设我们观察到一个结果为正面的硬币，现在它不再如之前所述的完全是红色(硬币B产生)和绿色(硬币C产生))，而是有一部分红色(概率是\(\mu_i\))，一部分绿色(概率是\(1-\mu_i\))，如下图所示

接下来我们就可以根据\(\mu_i\)重新估计参数\(\theta=(\pi,p,q)\)了，注意，为什么说是重新估计参数，因为\(\mu_i\)是根据当前给定参数\(\theta=(\pi,p,q)\)计算出来的：

根据之前我们估计参数的方法，有

\[\pi = \cfrac{红色的个数}{总个数} = \cfrac{\sum_{i=1}^N\mu_i}{N} \]

\[p=\cfrac{红色中1的个数}{红色的个数} =\cfrac{\sum_{i=1}^N\mu_ix_i}{\sum_{i=1}^N\mu_i} \]

\[q=\cfrac{绿色中1的个数}{绿色的个数}=\cfrac{\sum_{i=1}^N(1-\mu_i)x_i}{\sum_{i=1}^N(1-\mu_i)} \]

好了，通过不断的重复上述步骤，直到参数\(\theta=(\pi,p,q)\)不再变化，我们就得到了参数的最优估计。至此，我们直观地解释了EM算法。

E步：使用当前给定的\(\theta\)计算\(\mu\)，初始时可以根据先验指定\(\theta\)；

M步：使用E步计算得到的\(\mu\)重新计算\(\theta\)。

重复上述步骤，直至收敛。

3. EM算法推导

给定训练样本是\(\{x^{(1)},x^{(2)},\cdots,x^{(m)}\}\)，样例间独立，我们想要找到每个样例隐含类别\(z\)，能使得\(p(x,z)\)最大。\(p(x,z)\)的最大似然估计如下：

\[\begin{aligned} \ell(\theta) &= \sum_{i=1}^mlog \, p(x;\theta) \\ &= \sum_{i=1}^mlog \sum_{z}p(x,z;\theta) \end{aligned} \tag{3.1}\]

第一步是对极大似然取对数；第二步是对每个样例的每个可能类别\(z\)求联合分布概率和。但是直接求\(\theta\)一般比较困难，因为有隐藏变量\(z\)的存在，但是一般确定了\(z\)后，求解就容易了。

EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化\(\ell(\theta)\)，我们可以不断地建立\(\ell(\theta)\)的下界（E步），然后优化下界（M步）。这句话比较抽象，看下面的。

对于每个样例\(i\)，让\(Q_i\)表示该样例隐含变量\(z\)的某种分布，\(Q_i\)满足的条件是\(\sum_{z}Q_i(z)=1\)，\(Q_i(z) \ge 0\)。（如果\(z\)是连续性的，那么\(Q_i(z)\)是概率密度函数，需要将求和符号换做积分符号）。比如要将班上的学生聚类，假设隐藏变量\(z\)是身高，那么就是连续的高斯分布。如果按照隐藏变量\(z\)是性别，那就是伯努利分布了。

可以由前面阐述的内容得到下面的公式：

\[\begin{aligned} \sum_{i}log\, p(x^{(i)};\theta) &= \sum_{i}log \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta) \\ &= \sum_{i}log \sum_{z^{(i)}} Q_i(z^{(i)}) \cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\ &\ge \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \end{aligned} \tag{3.2}\]

此处推导利用Jensen不等式，考虑\(log(x)\)是凹函数，且

\[\sum_{z^{(i)}} Q_i(z^{(i)}) \bigg[\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\bigg] \]

可以看做\(\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\)的期望，可以看一下前面回顾的期望的公式。

对应于上述问题，根据式\((1.1)\)，\(Y\)是\(\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\)，\(X\)是\(z^{(i)}\)，\(Q_i(z^{(i)})\)是\(p_k\)，\(g\)是\(z^{(i)}\)到\(\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\)的映射。这样解释了期望的描述，再根据凹函数的Jensen不等式：

\[f\bigg(E_{z^{(i)\sim Q_i}}\bigg[\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\bigg]\bigg) \ge E_{z^{(i)\sim Q_i}}\bigg[f\bigg(\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\bigg)\bigg] \]

此时，\(f\)即是\(log\)函数，于是推导得到式\((3.2)\)。这个过程可以看做是对\(\ell(\theta)\)求下界函数。

现在，我们将式\((3.2)\)写成：对数似然函数\(\ell(\theta) \ge J(Q,\theta)\)，显然\(J(Q,\theta)\)是容易通过求导数，然后求得极值。但是由于是不等式，所以\(J(Q,\theta)\)的极大值并不是\(\ell(\theta)\)极大值，而我们想得到的是\(\ell(\theta)\)的极大值，那怎么办呢？

现在我们就需要一点想象力了，我们可以通过不断的最大化这个下界\(J\)，来使得\(\ell(\theta)\)不断提高，最终达到它的最大值。如下图所示

对于\(Q_i\)的选择，有多种可能，那么哪种更好呢？假设\(\theta\)已经给定，即\(\theta\)已经固定，那么\(\ell(\theta)\)的值就决定于\(Q_i(z^{(i)})\)和\(p(x^{(i)},z^{(i)})\)了。我们可以通过调整这两个概率使下界\(J(Q,\theta)\)不断上升，以逼近\(\ell(\theta)\)的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率值能够等价于\(\ell(\theta)\)了（浅蓝色曲线到蓝色曲线）。然后固定\(Q(z)\)，调整\(\theta\)使下界\(J(Q,\theta)\)达到最大值（\(\theta^t\)到\(\theta^{t+1}\)），重复上述过程直至收敛到对数似然函数\(\ell(\theta)\)的极大值处\(\theta^ \ast\)。

按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想等式成立，需要让随机变量变成常数值，这里得到：

\[\cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} = c \tag{3.3} \]

其中，\(c\)为常数，不依赖于\(z^{(i)}\)。对此式子做进一步推导，我们知道\(\sum_{z}Q_i(z^{(i)})=1\)，那么也就有\(\sum_{z}p(x^{(i)},z^{(i)};\theta)=c\)，（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是\(c\)），那么有下式：

\[\begin{aligned} Q_i(z^{(i)}) &= \cfrac{p(x^{(i)},z^{(i)};\theta)}{\sum_{z}p(x^{(i)},z^{(i)};\theta)} \\ & = \cfrac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &= p(z^{(i)}|x^{(i)};\theta) \end{aligned} \tag{3.4}\]

至此，我们推出了在固定其他参数\(\theta\)后，\(Q_i(z^{(i)})\)的计算公式就是后验概率，解决了\(Q_i(z^{(i)})\)如何选择的问题。这一步就是\(E\)步，建立\(\ell(\theta)\)的下界。接下来的\(M\)步，就是在给定\(Q_i(z^{(i)})\)后，调整\(\theta\)，极大化\(\ell(\theta)\)的下界（在固定\(Q_i(z^{(i)})\)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

循环重复直至算法收敛：

（\(E\)步）对于每一个样例\(i\)，计算

\[Q_i(z^{(i)}) := p(z^{(i)}|x^{(i)};\theta) \tag{3.5} \]

（\(M\)步）计算

\[\theta := arg \smash{\max_\theta} \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \tag{3.6} \]

那么，究竟怎么确保EM收敛？假定\(\theta^{(t)}\)和\(\theta^{(t+1)}\)是EM的第\(t\)次和\(t+1\)次迭代后的结果。如果我们证明了\(\ell(\theta^{(t)}) \le \ell(\theta^{(t+1)})\)，也就是说极大似然估计单调增加，那么最终我们会得到极大似然估计的最大值。下面来证明，选定\(\theta^{(t)}\)后，我们得到E步

\[Q_i^{(t)}(z^{(i)}) := p(z^{(i)}|x^{(i)};\theta^{(t)}) \]

这一步保证了在给定\(\theta^{(t)}\)时，Jensen不等式中的等式成立，也就是

\[\ell(\theta^{(t)})=\sum_{i} \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \]

然后进行M步，固定\(Q_i^{(t)}(z^{(i)})\)，并将\(\theta^{(t)}\)视作变量，对上面的\(\theta^{(t)}\)求导后，得到\(\theta^{(t+1)}\)，这样经过一些推导会有以下式子成立：

\[\begin{aligned} \ell{(\theta^{(t)})} &\ge \sum_{i} \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \\ &\ge \sum_{i} \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\ &= \ell(\theta^{(t)}) \end{aligned}\]

解释上式第一步，得到\(\theta^{(t+1)}\)时，只有最大化\(\ell(\theta^{(t)})\)，也就是\(\ell(\theta^{(t+1)})\)的下界，而没有使等式成立，等式成立只有在固定\(\theta\)，并按\(E\)步得到\(Q_i\)时才能成立。

况且根据我们前面得到的下式，对于所有的\(Q_i\)和\(\theta\)都成立

\[\ell(\theta) \ge \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \]

上式第二步，利用\(M\)步的定义，\(M\)步就是将\(\theta^{(t)}\)调整到\(\theta^{(t+1)}\)，使得下界最大化。

这样就证明了\(\ell(\theta\)会单调增加。一种收敛方法是\(\ell(\theta)\)不再变化，还有一种就是变化幅度很小。

再次解释一下。首先第一步对所有的参数都满足，而其等式成立条件只是在固定\(\theta\)，并调整好\(Q\)时成立，而第一步只是固定\(Q\)，调整\(\theta\)，不能保证等式一定成立。第一步到第二步就是M步的定义，第二步到第三步是前面E步所保证等式成立条件。也就是说E步会将下界拉到与\(\ell(\theta)\)一个特定值（这里\(\theta^{(t)}\)）一样的高度，而此时发现下界仍然可以上升，因此经过M步后，下界又被拉升，但达不到与\(\theta^{(t)}\)另外一个特定值一样的高度，之后E步又将下界拉到与这个特定值一样的高度，重复下去，直到最大值。

如果我们定义

\[J(Q,\theta) = \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) log \cfrac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \tag{3.7} \]

从前面的推导中我们知道\(\ell(\theta) \ge J(Q,\theta)\)，EM可以看做是\(J\)的坐标上升法，E步固定\(\theta\)优化\(Q\)，M步固定\(Q\)优化\(\theta\)。

4. 混合高斯模型(Gaussian mixture model, GMM)

4.1 单高斯模型(Gaussian single model, GSM）

简单回顾一下概率论讲过的高斯模型。

高斯模型是一种常用的随机变量分布模型，一维高斯分布的概率密度函数如下：

\[f(x) = \cfrac{1}{\sqrt{2\pi}\sigma}exp(-\cfrac{(x-\mu)^2}{2\sigma^2}) \tag{4.1} \]

\(\mu\)和\(\sigma^2\)是高斯分布的参数，分别表示均值和方差。

\(n\)维随机变量\(x=(x_1,x_2,\cdots,x_n)^T\)的联合概率密度函数为：

\[f(x) = \cfrac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}}exp[-\cfrac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)] \tag{4.2} \]

其中，\(n\)维向量\(\mu = (\mu_1,\mu_2,\cdots,\mu_n)^T\)被称为均值，\(n \times n\)的矩阵\(\Sigma\)被称为协方差，\(|\Sigma|\)表示\(\Sigma\)的行列式。

4.2 混合高斯模型

4.2.1 为什么要有混合高斯模型

先来看一组数据，如下图。

图1 某分布所产生数据

如果我们假设这组数据是由某个单一的高斯分布产生的，利用极大似然估计（后文会提到）对这个高斯分布做参数估计，得到一个最佳的高斯分布模型如下：

图2 用单高斯模型对样本作分析示意图

有什么问题吗？一般来讲越靠近椭圆的中心样本出现的概率越大，这是由概率密度函数决定的，但是这个高斯分布的椭圆中心的样本量却极少。显然样本服从单高斯分布的假设并不合理。单高斯模型无法产生这样的样本。

实际上，这是用两个不同的高斯分布混合产生的数据

图3 混合高斯模型对样本作分析示意图

通过求解两个高斯模型，并通过一定的权重将两个高斯模型融合成一个模型，即最终的混合高斯模型。这个混合高斯模型可以产生这样的样本。

更一般化的描述为：假设混合高斯模型有\(K\)个高斯模型组成，则GMM的概率密度函数如下：

\[p(x) = \sum_{k=1}^K p(k)p(x|k) = \sum_{k=1}^K \pi_k \mathcal{N}(x ;\mu_k,\Sigma_k) \tag{4.3} \]

其中，\(p(x|k)=\mathcal{N}(x ;\mu_k,\Sigma_k)\)是第\(k\)个高斯模型的概率密度函数，可以看成以\(k\)为条件\(x\)的概率；\(p(k) = \pi_k\)是第\(k\)个高斯模型的权重或称为混合系数，我们把其看成选择第\(k\)个模型的先验概率，且满足\(\sum_{k=1}^K\pi_k =1\)，\(0 \le \pi_k \le 1\)。后验概率\(p(k|x)\)起着一个重要的作用，它也被称为响应度(responsibilities)，表示对于每个样本\(x\)，它由第\(k\)个组份生成的概率，根据贝叶斯定理，后验概率可以表示为：

\[\begin{aligned} \gamma_k(x) &\equiv p(k|x) \\ &= \cfrac{p(k)p(x|k)}{\sum_{l=1}^K p(l)p(x|l)} \\ &= \cfrac{\pi_k \mathcal{N}(x ;\mu_k,\Sigma_k)}{\sum_{l=1}^K \pi_l \mathcal{N}(x ;\mu_l,\Sigma_l)} \end{aligned} \tag{4.4}\]

所以，混合高斯模型并不是什么新奇的东西，它的本质就是融合几个单高斯模型，来使得模型更加复杂，从而产生更复杂的样本。理论上，如果某个混合高斯模型融合的高斯模型个数足够多，它们之间的权重设定得足够合理，这个混合模型可以拟合任意分布的样本。

4.2.2 直观上理解混合高斯模型

下面通过几张图片来帮助理解混合高斯模型。

首先从简单的一维混合高斯模型说起。

图4 一维混合高斯模型示意图

在图4中，\(y1\),\(y2\)和\(y3\)分别表示三个一维高斯模型，他们的参数设定如图所示。\(y4\)表示将三个模型的概率密度函数直接相加，需要注意的是这并不是一个混合高斯模型，因为不满足\(\sum_{k=1}^K\pi_k =1\)的条件。而\(y5\)和\(y6\)分别是由三个相同的高斯模型融合生成的不同混合模型。由此可见，调整权重将极大影响混合模型的概率密度函数曲线。另一方面也可以直观地理解混合高斯模型可以更好地拟合样本的原因：它有更复杂更多变的概率密度函数曲线。理论上，混合高斯模型的概率密度函数曲线可以是任意形状的非线性函数。

下面再给出一个二维空间3个高斯模型混合的例子。

图5a 3个类别高斯分布截面轮廓线示意图

图5b 混合高斯分布截面轮廓线示意图

图5c 二维混合高斯分布概率密度函数示意图

图5a表示的是3个高斯模型的截面轮廓图，3个模型的权重系数已在图中注明，由截面轮廓图可知3个模型之间存在很多重叠区域。其实这也正是混合高斯模型所希望的。因为如果它们之间的重叠区域较少，那么生成的混合高斯模型一般较为简单，难以生成较为复杂的样本。

设定好了3个高斯模型和它们之间的权重系数之后，就可以确定二维混合高斯分布的概率密度函数曲面，如图5c所示。图5b是对于图5c概率密度曲面的截面轮廓线。从图5c也可以看出，整个混合高斯分布曲面相对比于单高斯分布曲面已经异常复杂。实际上，通过调整混合高斯分布的系数\(\pi, \mu, \Sigma\)，可以使得图5c的概率密度曲面去拟合任意的三维曲面，从而采样生成所需要的数据样本。

5. 混合高斯模型的求解

5.1 高斯模型的极大似然估计

假设我们采样得到随机变量\(X\)的一组样本\(\{x_1,x_2,\cdots,x_N\}\)，且我们认为变量\(X\)服从高斯分布，为了后续书写方便这里假设样本是一维的。数学形式表示为\(X \sim \mathcal{N}(\mu, \sigma)\)。采样的样本如图6所示，我们的目的就是找到一个合适的高斯分布，也就是确定高斯分布的参数\(\mu\)和\(\Sigma\)，使得这个高斯分布能产生这组样本的可能性尽可能大。

图6 最大化似然函数的意义

那怎么找到这个合适的高斯分布呢（在图8中的表示就是1~4哪个分布较为合适）？这时候似然函数就闪亮登场了。

下面我们就将上述问题数学化，并通过极大似然估计求解高斯模型的参数。令\(p(x_i;\mu,\sigma)\)是高斯分布的概率分布函数，表示\(X=x_i\)的概率。假设样本的抽样是独立的，那么我们同时抽到这\(N\)个样本的概率是抽到每个样本概率的乘积，也就是样本集\(X\)的联合概率。此联合概率即为似然函数：

\[L(\mu, \sigma) = \prod_{i=1}^N p(x_i;\mu,\sigma) = \prod_{i=1}^N \cfrac{1}{\sqrt{2\pi}\sigma}exp(-\cfrac{(x_i-\mu)^2}{2\sigma^2}) \]

先求对数似然函数，有

\[\ell(\mu,\sigma) = ln \, L(\mu, \Sigma) = \sum_{i=1}^N ln\bigg(\cfrac{1}{\sqrt{2\pi}\sigma}exp(-\cfrac{(x_i-\mu)^2}{2\sigma^2})\bigg) \]

接着对参数\(\mu, \sigma^2\)分别求导，并令其等于0

\[\cfrac{d\ell}{d\mu} = \sum_{i=1}^N\cfrac{(x_i-\mu)}{\sigma^2} = 0 \]

\[\cfrac{d\ell}{d\sigma^2} = \sum_{i=1}^N -\cfrac{1}{2\sigma^2} + \cfrac{(x_i-\mu)^2}{2\sigma^4}= 0 \]

解得

\[\mu = \cfrac{1}{N}\sum_{i=1}^Nx_i \qquad \sigma^2 = \cfrac{1}{N}\sum_{i=1}^N(x_i-\mu)^2 \]

上述结论所求出的参数就是最佳的高斯分布对应的参数，和矩估计的结果是一致的，并且意义非常直观：样本的均值即高斯分布的均值，样本的方差即高斯分布的方差。

所以最大化似然函数的意义就是：通过使得样本集的联合概率最大来对参数进行估计，从而选择最佳的分布模型。

对于图6产生的样本用极大化似然函数的方法，最终可以得到序号1对应的高斯分布模型是最佳的模型。

5.2 尝试用极大似然估计的方法来解GMM模型

高斯混合分布的形式由参数\(\pi, \mu, \Sigma\)控制。解GMM模型，实际上就是求得这组参数，使得其确定的GMM模型最有可能产生采样的样本。看上去和上一节的问题类似，我们我可以尝试利用极大似然估计的方法来求解参数。

如前所述，要利用极大似然估计求解模型最重要的一步就是求出似然函数，即样本集出现的联合概率。

根据混合高斯分布的公式，我们有对数似然函数

\[\ell(\pi,\mu,\Sigma) = ln\,p(X;\pi,\mu,\Sigma) = \sum_{i=1}^Nln\bigg(\sum_{k=1}^K\pi_k\mathcal{N(x_i;\mu_k,\Sigma_k)}\bigg) \tag{5.1} \]

好了，然后求导，令导数为0，得到模型参数\(\pi, \mu, \Sigma\)。

貌似问题已经解决了!

然而仔细观察可以发现，对数似然函数里面，对数里面还有求和。实际上没有办法通过求导的方法来求这个对数似然函数的最大值。

怎么办？

这时候就要要到EM算法求解。

如果我们已经清楚了某个变量服从的高斯分布，而且通过采样得到了这个变量的样本数据，想求高斯分布的参数，这时候极大似然估计可以胜任这个任务；而如果我们要求解的是一个混合模型，只知道混合模型中各个类的分布模型（譬如都是高斯分布）和对应的采样数据，而不知道这些采样数据分别来源于哪一类（隐变量），那这时候就可以借鉴EM算法。EM算法可以用于解决数据缺失的参数估计问题（隐变量的存在实际上就是数据缺失问题，缺失了各个样本来源于哪一类的记录）。

5.3 EM算法求解GMM

我们已经知道了EM的推导过程，这里就直接利用它解GMM。

E步：

按照一般EM公式得到

\[w_j^{(i)} = Q_i(z^{(i)}=j) = P(z^{(i)}=j|x^{(i)};\pi,\mu,\Sigma) \tag{5.2} \]

仔细观察，这不就是在描述混合高斯分布时说到的很重要的后验概率么，表示样本\(x_i\)由第\(j\)个组份生成的概率。

M步：

求得\(Q_i\)后之后，我们就得到了最贴近原对数似然函数的下界函数，那我们对它求极值就可以得到更新后的参数，先看一下这个下界函数

\[\begin{array}{l} \begin{aligned} \sum_{i=1}^N \sum_{z^{(i)}} Q_i(z^{(i)}) ln \cfrac{p(x^{(i)},z^{(i)};\pi, \mu, \Sigma)}{Q_i(z^{(i)})} \end{aligned}\\ \begin{aligned} &= \sum_{i=1}^N \sum_{j=1}^K Q_i(z^{(i)}=j) ln \cfrac{p(x^{(i)}|z^{(i)}=j;\mu, \Sigma)p(z^{(i)}=j;\pi)}{Q_i(z^{(i)}=j)} \\ &= \sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} ln \cfrac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp[-\cfrac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)]\cdot \pi_j}{w_j^{(i)}} \end{aligned} \end{array} \tag{5.3}\]

这是将\(z^{(i)}\)的\(k\)种情况展开后的样子，未知参数\(\pi_j\)，\(\mu_j\)和\(\Sigma_j\)。

固定\(\pi_j\)和\(\Sigma_j\)，对\(\mu_j\)求导得

\[\begin{array}{l} \begin{aligned} \nabla_{\mu_l}\sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} ln \cfrac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp[-\cfrac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)]\cdot \pi_j}{w_j^{(i)}} \end{aligned} \\ \begin{aligned} &= -\nabla_{\mu_l}\sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} \cfrac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j) \\ &= \cfrac{1}{2} \sum_{i=1}^N w_l^{(i)} \nabla_{\mu_l} 2 \mu_l^T\Sigma_l^{-1}x^{(i)} - \mu_l^T\Sigma_l^{-1}\mu_l \\ & = \sum_{i=1}^N w_l^{(i)}(\Sigma_l^{-1}x^{(i)}-\Sigma_l^{-1}\mu_l) \end{aligned} \end{array}\]

令上式等于0，得到

\[\mu_l = \cfrac{\sum_{i=1}^Nw_l^{(i)}x^{(i)}}{\sum_{i=1}^Nw_l^{(i)}} \tag{5.4} \]

这就是我们之前模型中的\(\mu\)的更新公式。

同样地，我们可以求出\(\Sigma\)的更新公式：

\[\Sigma_l = \cfrac{\sum_{i=1}^Nw_l^{(i)}(x^{(i)}-\mu_l)^T(x^{(i)}-\mu_l)}{\sum_{i=1}^Nw_l^{(i)}} \tag{5.5} \]

对于\(\pi\)的求解可能复杂一些，首先我们把下界函数中与\(\pi\)不相关的项全部去掉，实际上需要优化的公式是：

\[\sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} ln(\pi_j) \]

然后\(\pi\)作为各组份的混合系数有一个天然的约束条件就是所有的\(\pi_j\)之和为1，即

\[\sum_{j=1}^K \pi_j =1 \]

这个优化问题我们很熟悉了，直接构造拉格朗日乘子。

\[\mathcal{L}(\pi) = \sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} ln(\pi_j) + \beta(\sum_{j=1}^K \pi_j - 1) \]

还有一点就是\(\pi_j \ge 0\)，但这一点会在得到的公式里自动满足。求导得

\[\cfrac{\partial}{\partial \pi_j}\mathcal{L}(\pi) = \sum_{i=1}^N \cfrac{w_j^{(i)}}{\pi_j} + \beta \]

令其等于0，则

\[\pi_j = \cfrac{\sum_{i=1}^Nw_j^{(i)}}{-\beta} \]

再次使用\(\sum_{j=1}^K \pi_j =1\)，即将上式从1到\(K\)进行连加，有

\[-\beta = \sum_{i=1}^N \sum_{j=1}^K w_j^{(i)} = \sum_{i=1}^N 1 = N \]

这样就神奇地得到了\(\beta\)。

于是得到M步中\(\pi\)的更新公式：

\[\pi_j = \cfrac{1}{N}\sum_{i=1}^N w_j^{(i)} \tag{5.6} \]

我们再来看看这些参数的意义，其实未尝不符合我们的直觉认识。\(w_j^{(i)}\)可以看做第\(i\)个样本属于第\(j\)类的概率，那么所有样本中属于第\(j\)类的个数就是\(w_j^{(i)}\)之和，每个样本\(x_i\)对应第\(j\)类的值就是\(w_j^{(i)}x_i\)，这样算的平均数就是第\(j\)类对应的\(\mu\)，继续按照这个思路计算的方差就是第\(j\)类的\(\Sigma\)，第\(j\)类的概率就是第\(j\)类的个数除以总样本数。所以，GMM模型虽然推导起来有点吓人，但仔细想想它最后的结果也是符合我们的直觉认识的，每个样本都是一部分属于某一类，所有样本中的某一类的部分构成了这一类的分布。