EM算法之高斯混合模型（一）

单个高斯模型

如果我们有一堆数据，其分布属于一个高斯模型，那么有

p(X)=N(x|μ,Σ)=1(2π)m|Σ|‾‾‾‾‾‾‾‾√exp[−12(x−μ)TΣ−1(x−μ)](1.1)

这样子的话，对于单个高斯，我们可以直接对其参数μ和Σ进行求导，求出对应的参数。
那么现在有一堆数据，其分布如下所示，

那么我们需要用多个高斯对数据的分布进行描述。接下来我们看看多个高斯混合模型.

混合高斯模型

每个GMM由K个Gaussian分布组成，每个高斯分布(Gaussian)称为一个“Component”，这些Component 线性加成在一起就组成了 GMM 的概率密度函数：

p(x)=∑k=1Kp(k)p(x|k)=∑k=1KπkNk(x|μk,Σk)(2.1)

上式中∑Kk=1πk=1，其中：

Nk(x|μk,Σk)=1(2π)m|Σk|‾‾‾‾‾‾‾‾‾√exp[−12(x−μk)TΣ−1k(x−μk)](2.2)

这边我们的参数可以用θ⎯⎯来表示：
θ⎯⎯={μ1,μ2...μk,Σ1,Σ2,...,Σk,π1,π2,...,πk}

这边的πi表示的是每个高斯分布对数据分布的权重，∑kiπi=1

那我们如何从分布中取一点吗？我们可以分成两步，首先πk的概率选择一个component(每个component就是一个聚类中心)，然后再从选中的这个歌高斯分布中抽取一个点。

高斯混合模型参数估计与似然函数

如果我们直接按照单个高斯分布那样直接对高斯混合模型使用Maximum Likelihood来求解的话：

θ⎯⎯MLE=argmaxθ{∑iNln[∑ikπlN(μl,Σl)]}(3.1)

分别对每个高斯进行求导，式子中log里面的都是求和，这样子分别求导是很困难的。

EM算法进行参数求解

最近看到一篇深入讲解EM算法的文章（链接在最后），然后赶紧把其中的东西拿过来补充一下EM算法。
EM算法就是E 期望 + M 最大化两步。那么我们先看一个直观的例子：

最经典的例子就是抛3个硬币，跑I硬币决定C1和C2，然后抛C1或者C2决定正反面， 然后估算3个硬币的正反面概率值。

这个例子为什么经典， 因为它告诉我们，当存在隐变量I的时候， 直接的最大似然估计无法直接搞定。

EM算法
输入：观测数据X，隐变量数据Z，联合分布P(X,Z|θ),条件分布P(Z|X,θ);
输出模型参数：θ
1. 选择初始参数θ(0),开始迭代;
2. E step:记θ(i)为第i次迭代时θ的参数估计，那么第i+1步迭代记做:
Q(θ,θ(i))=Ez[log(X,Z|θ)|X,θ(i)]=∑zlogP(X,Z|θ)P(Z|X,θ(i))
这里的Q(θ,θ(i))是对数似然函数logP(X,Z|θ)关于在给定观测数据X和当前参数θ(i)下对未观测数据Z的条件概率分布P(Z|Y,θ(i))。这边如果隐含变量Z是连续的话，我们可以使用积分∫z而不是求和∑z
3. M step:求Q(θ,θ(i))的最大化,即为我们第i+1次迭代参数θ(i)的值。
θ(i+1)=argmaxθQ(θ,θ(i))
4. 重复第2和第3步，直到收敛。迭代停止的条件是，一般满足：
||θ(i+1)−θ(i)||<ε 或者 ||Q(θ(i+1),θ(i)||<ε

根据上述EM算法我们给出我们参数求解的迭代公式：

θ(g+1)=argmaxθ∫zln{P(x,z|θ)P(z|x,θ(g))}dz(4.1)

x : 是我们已有的数据
z : 是隐变量，如果隐变量是连续的就用积分，如果是离散的就用求和。

加入隐变量不能改变原有的边缘分布，即：
p(x) = ∫zP(x|z)P(z)dz

那这边的隐变量是什么呢？那我们先看一个图：

我们有一些数据其分布属于两个混合高斯，每个数据xi有一个对应的zi，这边的z就是我们的隐变量，而每个zi等于1或者2，1指的是该数据属于第一个高斯分布，2指的是该数据属于第二个高斯分布。
那我们看看加入隐变量会不会改变原有的分布：

P(xi)=∑ziPθ(xi|zi)Pθ(zi)=∑ziπziN(μzi,Σzi)(4.2)

Pθ(zi)=πzi 所以我们可以从上式中可以看出加入隐变量并未改变原有的分布。

EM算法的收敛性

设P(X|θ)为观测数据的似然函数，θ(i)为EM算法得到的参数估计序列，P(X|θ(i))为对应的似然函数序列，则P(X|θ(i))是单调递增的，即

P(X|θ(i+1))≥P(X|θ(i))(5.1)

proof:
要证明函数收敛，就是证明：

logP(X|θ(g+1))≥logP(X|θ(g))

那么

logP(X|θ)=logP(X,Z|θ)−log(Z|X,θ)(5.2)

因为P(X|θ)=P(X,Z|θ)P(Z|Y,θ)，即Pθ(X)=Pθ(X,Z)Pθ(Z|X)

现在我们对上面的式子分别对P(Z|Y,θ(i))求期望：

左边=EP(Z|X,θ(i)){logP(X|θ)}=∫zlog(P(X|θ))P(Z|X,θ(i))dz=log(P(X|θ))∫zP(Z|X,θ(i))=log(P(X|θ))(5.3)

那么右边

右边=∫zlogP(X,Z|θ)P(Z|X,θ(i))dz−∫zlogP(Z|X,θ)P(Z|X,θ(i))dz(5.4)

现在我们令

Q(θ,θ(i))=∫zlogP(X,Z|θ)P(Z|X,θ(i))dz=∑zlogP(X,Z|θ)P(Z|X,θ(i))(5.5)H(θ,θ(i))=∫zlogP(Z|X,θ)P(Z|X,θ(i))dz=∑zlogP(Z|X,θ)P(Z|X,θ(i))(5.6)

于是对数似然函数可以写成：

log(P(X|θ))=Q(θ,θ(i))−H(θ,θ(i))(5.7)

在式子(5.7)中分别取θ为θ(i)和θ(i+1):

logP(X|θ(i+1))−logP(X|θ(i))=[Q(θ(i+1),θ(i))−Q(θ(i),θ(i))]−[H(θ(i+1),θ(i))−H(θ(i),θ(i))](5.8)

因为根据(4.1)，θ(i+1)使Q(θ,θ(i))达到最大，所以有：

Q(θ(i+1),θ(i))−Q(θ(i),θ(i))≥0(5.9)

那

H(θ(i+1),θ(i))−H(θ(i),θ(i))=∑z⟮logP(Z|X,θ(i+1))P(Z|X,θ(i))⟯P(Z|X,θ(i))≤log⟮∑zP(Z|X,θ(i+1))P(Z|X,θ(i))⟯P(Z|X,θ(i))=log⟮P(Z|X,θ(i+1))⟯=0(5.9)

这边的不等式用到了jensen不等式，pf(x)+(1-p)f(y) ≥ f[(1-p)y + px]

所以综上EM算法收敛。

参考：

贝叶斯推断及其互联网应用（一）：定理简介

漫谈 Clustering (3): Gaussian Mixture Model

GMM算法（Python版）

EM算法的九层境界：​Hinton和Jordan理解的EM算法