机器学习——概率分类(三)高斯概率密度与混合高斯模型

在之前的文章机器学习——概率分类(一)朴素贝叶斯模型一文中，我们主要介绍了机器学习中的概率分类问题。我们提出了简单的朴素贝叶斯模型来进行概率密度的估计。在本篇文章中，我们主要介绍概率密度估计的第二种方法——高斯密度估计。

1 高斯概率密度原理

1.1 高斯概率密度引入

首先，我们假设样本X符合的是高斯分布，当X的维度是一维的时候，其符合的概率分布的公式为：
P ( X ) = 1 2 π δ e ( x − u ) 2 2 δ P(X)=\frac{1}{\sqrt{2πδ}}e^{\frac{(x-u)^2}{2δ}} P(X)=2πδ 1e2δ(x−u)2
u = 1 N ∑ i = 1 N X i u = \frac{1}{N}∑_{i=1}^NX_i u=N1i=1∑NXi
δ 2 = 1 N ∑ i = 1 N ( X i − u ) 2 δ^2=\frac{1}{N}∑_{i=1}^N(X_i-u)^2 δ2=N1i=1∑N(Xi−u)2
其中u是高斯分布的均值， δ 2 δ^2 δ2是高斯分布的方差。下面，当X是一个多维度的样本的时候，我们可以确定其高斯分布的公式为：
P ( X ) = 1 ( 2 π ) d ∣ ε ∣ e x p [ − 1 2 ( X − u ) T ε − 1 ( X − u ) ] P(X)=\frac{1}{\sqrt{(2π)^d|ε|}}exp[-\frac{1}{2}(X-u)^Tε^{-1}(X-u)] P(X)=(2π)d∣ε∣ 1exp[−21(X−u)Tε−1(X−u)]
其中，X是一个d维度的向量，u是d维的向量，表示均值。其中ε表示d*d的协方差矩阵，协方差矩阵具体可以表示为：
ε = ε 11 ε 12 . . . ε 1 d ε 21 ε 22 . . . ε 2 d … … … … ε d 1 ε d 2 . . . ε d d ε= \begin{matrix} ε_{11}&&ε_{12}&&...&&ε_{1d}\\ ε_{21}&&ε_{22}&&...&&ε_{2d}\\ &&…………\\ ε_{d1}&&ε_{d2}&&...&&ε_{dd}\\ \end{matrix} ε=ε11ε21εd1ε12ε22…………εd2.........ε1dε2dεdd
其中任意的 ε i j = E [ ( X i − E ( X i ) ) ( X j − E ( X j ) ) ] ε_{ij}=E[(X_i-E(X_i))(X_j-E(X_j))] εij=E[(Xi−E(Xi))(Xj−E(Xj))]。并且，ε矩阵是对称的。

1.2 高斯概率密度推导

根据上面的描述，我们知道了，存在一个样本集合D，其中包含N个样本X，这些样本符合高斯分布。那么我们现在要做的就是推导是该高斯分布的u和δ。进而我们就能够确定各个样本的概率密度。

下面我们利用极大似然算法来设计目标函数：
E ( u , ε ) = ∑ i = 1 N l n ( P ( X i ) ) E(u,ε)=∑_{i=1}^Nln(P(X_i)) E(u,ε)=i=1∑Nln(P(Xi))
其中E(u,ε)是待定的函数，u，ε是函数中待定的参数。 P ( X i ) P(X_i) P(Xi)是 X i X_i Xi的高斯密度函数，C表示类别。之所以加上ln，是因为高斯密度函数中exp的存在，通过ln可以简单运算。

现在，我们给出两个假设：

假设1：所有的样本X是独立同分布的。

假设2：求得的参数u和ε是的{ X i X_i Xi},i∈[1,N]的概率最大。

下一步，我们将高斯密度函数带入到我们设定的目标函数中，有：
E ( u , ε ) = ∑ i = 1 N l n ( 1 ( 2 π ) d ∣ ε ∣ e x p [ − 1 2 ( X − u ) T ε − 1 ( X − u ) ] ) E(u,ε)=∑_{i=1}^Nln(\frac{1}{\sqrt{(2π)^d|ε|}}exp[-\frac{1}{2}(X-u)^Tε^{-1}(X-u)]) E(u,ε)=i=1∑Nln((2π)d∣ε∣ 1exp[−21(X−u)Tε−1(X−u)])
化简之后得到：
E ( u , ε ) = − N d 2 l n ( 2 π ) − N 2 l n ( ∣ ε ∣ ) − 1 2 ∑ i = 1 N ( X i − u ) T ε − 1 ( X i − u ) E(u,ε)=-\frac{Nd}{2}ln(2π)-\frac{N}{2}ln(|ε|)-\frac{1}{2}∑_{i=1}^N(X_i-u)^Tε^{-1}(X_i-u) E(u,ε)=−2Ndln(2π)−2Nln(∣ε∣)−21i=1∑N(Xi−u)Tε−1(Xi−u)

既然，我们要求极值，那么自然而然的能够想到的就是求导：

∂ E ∂ u = − ∑ i = 1 N ε − 1 ( X i − u ) \frac{ ∂E}{ ∂u}=-∑_{i=1}^Nε^{-1}(X_i-u) ∂u∂E=−i=1∑Nε−1(Xi−u)
∂ E ∂ ε = − N 2 ε − 1 + 1 2 ε − 1 { ∑ i = 1 N ( X i − u ) ( X i − u ) T } ε − 1 \frac{ ∂E}{ ∂ε}=-\frac{N}{2}ε^{-1}+\frac{1}{2}ε^{-1}\{∑_{i=1}^N(X_i-u)(X_i-u)^T\}ε^{-1} ∂ε∂E=−2Nε−1+21ε−1{i=1∑N(Xi−u)(Xi−u)T}ε−1
这里没有给出公式具体的推导过程，具体可以参考多维高斯概率密度函数对协方差矩阵求导这篇文章。

下面就是最后一步了，我们让u和 ε ε ε导数为0，就可以求出极值。
∂ E ∂ u = − ∑ i = 1 N ε − 1 ( X i − u ) = 0 ，则 ∑ i = 1 N ( X i − u ) = 0 ，则 u = 1 N ∑ i = 1 N X i \frac{ ∂E}{ ∂u}=-∑_{i=1}^Nε^{-1}(X_i-u)=0，则∑_{i=1}^N(X_i-u)=0，则u=\frac{1}{N}∑_{i=1}^NX_i ∂u∂E=−i=1∑Nε−1(Xi−u)=0，则i=1∑N(Xi−u)=0，则u=N1i=1∑NXi
∂ E ∂ ε = − N 2 ε − 1 + 1 2 ε − 1 { ∑ i = 1 N ( X i − u ) ( X i − u ) T } ε − 1 = 0 \frac{ ∂E}{ ∂ε}=-\frac{N}{2}ε^{-1}+\frac{1}{2}ε^{-1}\{∑_{i=1}^N(X_i-u)(X_i-u)^T\}ε^{-1}=0 ∂ε∂E=−2Nε−1+21ε−1{i=1∑N(Xi−u)(Xi−u)T}ε−1=0
在等式的两侧分别乘以ε，则有
N 2 = 1 2 ε − 1 ∑ i = 1 N ( X i − u ) ( X i − u ) T \frac{N}{2}=\frac{1}{2}ε^{-1}∑_{i=1}^N(X_i-u)(X_i-u)^T 2N=21ε−1i=1∑N(Xi−u)(Xi−u)T
解出：
ε = 1 N ∑ i = 1 N ( X i − u ) ( X i − u ) T ε=\frac{1}{N}∑_{i=1}^N(X_i-u)(X_i-u)^T ε=N1i=1∑N(Xi−u)(Xi−u)T

总结：

通过求导过程，我们最终确定了该高斯分布的两个参数：
u = 1 N ∑ i = 1 N X i u=\frac{1}{N}∑_{i=1}^NX_i u=N1i=1∑NXi
ε = 1 N ∑ i = 1 N ( X i − u ) ( X i − u ) T ε=\frac{1}{N}∑_{i=1}^N(X_i-u)(X_i-u)^T ε=N1i=1∑N(Xi−u)(Xi−u)T

2 混合高斯模型GMM

2.1 GMM模型的引入

在我们之前的高斯概率密度估计中，我们设定的是所有的样本的都是属于同一个高斯分布的，我们的目标是确定这种高斯分布。类似于如下的情况：

上图所示的是一个二维的平面，我们可以看出的是，所有的数据都接近于一个高斯分布。但是，作为一个分类问题而言，数据的分布肯定是有区别的，否则同属于一个分布，就无法处理分类问题了。所以数据的分布肯定是有差异的，我们拿下一个图举一个例子：

我们可以看到的是，所以的样本数据可以大致本分成两类，也就是说，存在两个概率密度的分布(由红色和蓝色的圈表示)。显然我们不可能在用单一的高斯分布来计算所所有的数据的概率密度函数，那么我们应该如何处理呢，这就引出来了我们下面要介绍的模型——高斯混合模型。

2.2 GMM模型的原理推导

现在假设，我们的数据样本集合D一共包含K个类别。每个类别k下的样本均符合高斯分布，则我们在考虑不同样本的概率密度的时候，需要考虑它在不同分类下的高斯分布的概率密度，我们可以做一个加权求和来确定其最终的概率密度函数。数学化表示一下就是：
P ( X ) = ∑ k = 1 K Π k N ( X ∣ u k , ε k ) P(X)=∑_{k=1}^KΠ_kN(X|u_k,ε_k) P(X)=k=1∑KΠkN(X∣uk,εk)
N ( X ∣ u k , ε k ) = 1 ( 2 π ) d ∣ ε k ∣ e x p [ − 1 2 ( X − u k ) T ε k − 1 ( X − u k ) ] N(X|u_k,ε_k)=\frac{1}{\sqrt{(2π)^d|ε_k|}}exp[-\frac{1}{2}(X-u_k)^Tε_k^{-1}(X-u_k)] N(X∣uk,εk)=(2π)d∣εk∣ 1exp[−21(X−uk)Tεk−1(X−uk)]
其中N(X|u_k,ε_k)表示第k个分类的高斯密度，也就是第k个分类产生X的概率。 Π k Π_k Πk表示第k个分类

的权重，其满足 ∑ k = 1 K Π k = 1 ∑_{k=1}^KΠ_k=1 ∑k=1KΠk=1。这样，我们就可以根据不同的k类来表示不同样本的x的最终的

概率密度函数。

下面，我们继续利用极大似然算法来假设目标函数：
E ( { Π k , u k , ε k } ) = ∑ i = 1 N l n ( ∑ k = 1 K Π k N ( X ∣ u k , ε k ) ) E(\{Π_k,u_k,ε_k\})=∑_{i=1}^Nln(∑_{k=1}^KΠ_kN(X|u_k,ε_k)) E({Πk,uk,εk})=i=1∑Nln(k=1∑KΠkN(X∣uk,εk))
将 N ( X ∣ u k , ε k ) ) N(X|u_k,ε_k)) N(X∣uk,εk))带入到函数中后有：
E ( { Π k , u k , ε k } ) = ∑ i = 1 N l n ( ∑ k = 1 K Π k 1 ( 2 π ) d ∣ ε k ∣ e x p [ − 1 2 ( X − u k ) T ε k − 1 ( X − u k ) ] ) E(\{Π_k,u_k,ε_k\})=∑_{i=1}^Nln(∑_{k=1}^KΠ_k\frac{1}{\sqrt{(2π)^d|ε_k|}}exp[-\frac{1}{2}(X-u_k)^Tε_k^{-1}(X-u_k)]) E({Πk,uk,εk})=i=1∑Nln(k=1∑KΠk(2π)d∣εk∣ 1exp[−21(X−uk)Tεk−1(X−uk)])

我们的最终要求的是目标函数的极值，常规来讲，也就是要目标函数的导数为0。数学化表示就是：
∂ E ∂ u k = 0 ， ∂ E ∂ π k = 0 ， ∂ E ∂ ε k = 0 \frac{ ∂E}{ ∂u_k}=0，\frac{ ∂E}{ ∂π_k}=0，\frac{ ∂E}{ ∂ε_k}=0 ∂uk∂E=0，∂πk∂E=0，∂εk∂E=0
对于目标函数，由于其是一个非凸的函数，那么我们在求解问题的时候，就不可能之间利用导数。求取该目标函数的方法主要包括三种，梯度下降算法，启发式算法，EM算法。这里我们选择EM算法来求解目标。

2.3 EM算法求解高斯混合模型

2.3.1 为什么选择EM算法求解？

对于EM算法而言，其有以下几个特点：

EM算法求解的是局部极值，而非全局极值。
EM算法无需调参，其一定是收敛的。
EM算法的实现方式简单，理论公式优美。
EM算法只针对某一类的局部极值问题求解，如下图所示：

如上图所示的样本是一个一维分布，我们可以很明显的看出来，样本大致属于两个高斯分布，但是我们不能确定两个高斯分布的均值也方差。此时，我们可以利用EM算法来估计两个分布的均值与方差。

2.3.2 求解流程

随机给各个样本X指定一个分类。在这种情况下，如果X属于 c i c_i ci类别的时候，有 P ( X ∣ c i ) > P ( X ∣ c j ) P(X|c_i)>P(X|c_j) P(X∣ci)>P(X∣cj)。
随机初始化 u k ， π k ， ε k ， k ∈ [ 1 , K ] u_k，π_k，ε_k，k∈[1,K] uk，πk，εk，k∈[1,K]。
E-STEP：
求解：
γ i k = π k N ( X i ∣ u k , ε k ) ∑ k = 1 K π k N ( X i ∣ u k , ε k ) ， k ∈ [ 1 , K ] ， i ∈ [ 1 , N ] γ_{ik}=\frac{π_kN(X_i|u_k,ε_k)}{∑_{k=1}^Kπ_kN(X_i|u_k,ε_k)}，k∈[1,K]，i∈[1,N] γik=∑k=1KπkN(Xi∣uk,εk)πkN(Xi∣uk,εk)，k∈[1,K]，i∈[1,N]
也就是说，我们分别计算每一个样本 X i X_i Xi在不同分类下的高斯分布概率与其总的高斯分布概率密度的比值。
M-STEP:
在这一步中，为了简化运算，我们首先定义两个基本的公式:
Q k = ∑ i = 1 N γ i k ， Q = ∑ k = 1 K Q k Q_k=∑_{i=1}^Nγ_{ik}，Q=∑_{k=1}^KQ_k Qk=i=1∑Nγik，Q=k=1∑KQk
第一个公式说明的是所有样本在类别k下面的高斯密度比值之和。第二个公式是将所有的类别的比值和再次累加，得到的是所有样本关于所有分类的高斯权重比值之和。有了上面的两个公式，我们就可以来更新相关的参数值了：
π k n e w = N k N π_k^{new}=\frac{N_k}{N} πknew=NNk
u k n e w = 1 N k ∑ i = 1 N γ i k ∗ X i u_k^{new}=\frac{1}{N_k}∑_{i=1}^Nγ_{ik}*X_i uknew=Nk1i=1∑Nγik∗Xi
ε k n e w = 1 N k ∑ i = 1 N γ i k ( X i − u k n e w ) ( X i − u k n e w ) T ε_k^{new}=\frac{1}{N_k}∑_{i=1}^Nγ_{ik}(X_i-u_k^{new})(X_i-u_k^{new})^T εknew=Nk1i=1∑Nγik(Xi−uknew)(Xi−uknew)T
从第4步回到第3步继续计算，知道算法收敛结束。