EM算法在高斯混合模型学习中的应用

本篇文章是之前期望极大算法(EM算法)文章的后续，有需要可以先看看那篇文章关于EM算法的推导。

高斯混合模型

高斯混合模型是研究算法的人避不开的一个东西，其在非深度学习的远古时代经常被用到，比如图像处理任务的前背景提取，点云处理任务的点云聚类等等等等。

具体的，高斯混合模型是指具有如下形式的概率分布模型：
P(y∣θ)=∑k=1Kαkϕ(y∣θk)P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right) P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中，αk\alpha_{k}αk是系数，αk⩾0\alpha_{k} \geqslant 0αk⩾0,∑k=1Kαk=1\quad \sum_{k=1}^{K} \alpha_{k}=1∑k=1Kαk=1；ϕ(y∣θk)\phi\left(y \mid \theta_{k}\right)ϕ(y∣θk)是高斯分布密度，θk=(μk,σk2)\theta_{k}=\left(\mu_{k}, \sigma_{k}^{2}\right)θk=(μk,σk2)，
ϕ(y∣θk)=12πσkexp⁡(−(y−μk)22σk2)\phi\left(y \mid \theta_{k}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right) ϕ(y∣θk)=2πσk1exp(−2σk2(y−μk)2)
称为第kkk个分模型。

QQQ函数的一般表达

在算法处理的过程中，将问题建模成高斯混合模型后，往往需要去解模型中的参数，这个时候就需要用到EM算法。

在期望极大算法(EM算法)文章的分析中我们已经知道，要进行EM算法，得先得到QQQ函数：
Q(θ,θ(i))=∑Zlog⁡P(Y,Z∣θ)P(Z∣Y,θ(i))Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} \log P(Y, Z \mid \theta)P\left(Z \mid Y, \theta^{(i)}\right) Q(θ,θ(i))=Z∑logP(Y,Z∣θ)P(Z∣Y,θ(i))
在随后的抛硬币问题中我们也介绍了一种求解这个QQQ函数的方法.但如果看过李老师的书的人会发现，本人没有用书中给出的看起来更具总结性的QQQ函数形式来求解问题，原因在于当时本人也不明白那个式子是怎么来的。。。

在李航老师的书上，QQQ函数是这样定义的：完全数据的对数似然函数log⁡P(Y,Z∣θ)\log P(Y, Z \mid \theta)logP(Y,Z∣θ)关于在给定观测数据YYY和当前参数θ(i)\theta^{(i)}θ(i)下对未观测数据ZZZ的条件概率分布P(Z∣Y,θ(i))P\left(Z \mid Y, \theta^{(i)}\right)P(Z∣Y,θ(i))的期望称为QQQ函数，即：
Q(θ,θ(i))=EZ[log⁡P(Y,Z∣θ)∣Y,θ(i)]Q\left(\theta, \theta^{(i)}\right)=E_{Z}\left[\log P(Y, Z \mid \theta) \mid Y, \theta^{(i)}\right] Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]
第一次看到这个定义和下面的公式感觉整个人都不好了！实在不知道他们之间为什么是个相等的关系。在经过一两个小时的发呆、无助、掉头发后突然看懂了！

假设log⁡P(Y,Z∣θ)\log P(Y, Z \mid \theta)logP(Y,Z∣θ)是只与变量ZZZ相关的函数，则可以把其写成f(Z)f(Z)f(Z)，当ZZZ取得一个定值的时候，其就是一个固定的数值。如果这个时候对它取期望，就有：
EZ(f(Z))=∑Z[f(Z)P(Z)]E_Z(f(Z))=\sum_{Z}\left[f(Z)P\left(Z\right)\right] EZ(f(Z))=Z∑[f(Z)P(Z)]
而如果Z的取值本身受到别的参数xxx，yyy影响，而这些参数都已经给出，则原式可以写成：
E(f(Z)∣x,y)=∑Z[f(Z)P(Z∣x,y)]E(f(Z)\mid x,y)=\sum_{Z}\left[f(Z)P\left(Z\mid x,y\right )\right] E(f(Z)∣x,y)=Z∑[f(Z)P(Z∣x,y)]
代入我们已知的量，等式得证。

由证明的过程也可以知道，这里的给定观测数据YYY和当前参数两个已知量影响的只有ZZZ这一变量。由于他们与log⁡P(Y,Z∣θ)\log P(Y, Z \mid \theta)logP(Y,Z∣θ)中的两个参数看起来似乎有关系，因此才大大增加了理解的难度。

有了QQQ函数的一般表达，从式子的形式我们知道关键的一步是把高斯混合模型的log⁡P(Y,Z∣θ)\log P(Y, Z \mid \theta)logP(Y,Z∣θ)给列出来，也就是把其完全数据的似然函数的对数列出来。

高斯混合模型参数估计的EM算法

假设数据y1,y2,⋯,yNy_{1},y_{2}, \cdots, y_{N}y1,y2,⋯,yN由高斯混合模型生成，
P(y∣θ)=∑k=1Kαkϕ(y∣θk)P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right) P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中，θ=(α1,α2,⋯,αK;θ1,θ2,⋯,θK)\theta=\left(\alpha_{1}, \alpha_{2}, \cdots, \alpha_{K} ; \theta_{1}, \theta_{2}, \cdots, \theta_{K}\right)θ=(α1,α2,⋯,αK;θ1,θ2,⋯,θK)，求高斯混合模型的参数，就是用EM算法估计高斯混合模型的参数θ\thetaθ。

回顾观测数据yj,j=1,2,⋯,N,y_{j},j=1,2, \cdots, N,yj,j=1,2,⋯,N, 的产生过程：

依概率αk\alpha_{k}αk选择第kkk个高斯分布分模型ϕ(y∣θk)\phi\left(y \mid \theta_{k}\right)ϕ(y∣θk)；
依第kkk个分模型的概率分布ϕ(y∣θk)\phi\left(y \mid \theta_{k}\right)ϕ(y∣θk)生成观测数据yjy_{j}yj。

在高斯混合模型建模中，结果通常是已知的—观测数据yj,j=1,2,⋯,N,y_{j},j=1,2, \cdots, N,yj,j=1,2,⋯,N, 是已知的；

而结果产生自哪个分模型未知—反映观测数据yjy_{j}yj来自第kkk个分模型的数据未知，k=1,2,⋯,K，k=1,2, \cdots, K ，k=1,2,⋯,K，以隐变量γjk\gamma_{j k}γjk表示，可定义如下:
γjk={1,第 j个观测来自第 k个分模型 0,否则 \gamma_{j k}=\left\{\begin{array}{ll} 1, & \text { 第 } j \text { 个观测来自第 } k \text { 个分模型 } \\ 0, & \text { 否则 } \end{array}\right. γjk={1,0, 第 j 个观测来自第 k 个分模型否则

j=1,2,⋯,N;k=1,2,⋯,Kj=1,2, \cdots, N ; \quad k=1,2, \cdots, K j=1,2,⋯,N;k=1,2,⋯,K

有了观测数据yjy_jyj及未观测数据γjk\gamma_{j k}γjk，那么完全数据是
(yj,γj1,γj2,⋯,γjK),j=1,2,⋯,N\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}\right), \quad j=1,2, \cdots, N (yj,γj1,γj2,⋯,γjK),j=1,2,⋯,N
于是，可以写出完全数据的似然函数：
P(y,γ∣θ)=∏j=1NP(yj,γj1,γj2,⋯,γjK∣θ)=∏k=1K∏j=1N[αkϕ(yj∣θk)]γjk=∏k=1Kαkn∏j=1N[ϕ(yj∣θk)]γjk=∏k=1Kαknk∏j=1N[12πσkexp⁡(−(yj−μk)22σk2)]γjk\begin{array}{l} P(y, \gamma \mid \theta)=\prod_{j=1}^{N} P\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K} \mid \theta\right)\\ =\prod_{k=1}^{K} \prod_{j=1}^{N}\left[\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{jk}} \\ =\prod_{k=1}^{K} \alpha_{k}^{n} \prod_{j=1}^{N}\left[\phi\left(y_{j} \mid \theta_{k}\right)\right]^{\gamma_{jk}} \\ =\prod_{k=1}^{K} \alpha_{k}^{n_{k}} \prod_{j=1}^{N}\left[\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)\right]^{\gamma_{jk}} \end{array} P(y,γ∣θ)=∏j=1NP(yj,γj1,γj2,⋯,γjK∣θ)=∏k=1K∏j=1N[αkϕ(yj∣θk)]γjk=∏k=1Kαkn∏j=1N[ϕ(yj∣θk)]γjk=∏k=1Kαknk∏j=1N[2πσk1exp(−2σk2(yj−μk)2)]γjk
式中，nk=∑j=1Nγjkn_{k}=\sum_{j=1}^{N} \gamma_{j k}nk=∑j=1Nγjk，∑k=1Knk=N\sum_{k=1}^{K} n_{k}=N∑k=1Knk=N。

那么，完全数据的对数似然函数为：
log⁡P(y,γ∣θ)=∑k=1Knklog⁡αk+∑j=1Nγjk[log⁡(12πσk)−(yj−μk)22σk2]=∑k=1Knklog⁡αk+∑j=1Nγjk[log⁡(12π)−log⁡σk−(yj−μk)22σk2]\log P(y, \gamma \mid \theta)=\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}\sigma_{k}}\right)-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\\ =\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log{\sigma_{k}-}\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\\ logP(y,γ∣θ)=k=1∑Knklogαk+j=1∑Nγjk[log(2πσk1)−2σk2(yj−μk)2]=k=1∑Knklogαk+j=1∑Nγjk[log(2π1)−logσk−2σk2(yj−μk)2]
整个QQQ函数为：
Q(θ,θ(i))=E[log⁡P(y,γ∣θ)∣y,θ(i)]=E{∑k=1Knklog⁡αk+∑j=1Nγjk[log⁡(12π)−log⁡σk−(yj−μk)22σk2]}=∑k=1K{∑j=1N(Eγjk)log⁡αk+∑j=1N(Eγjk)[log⁡(12π)−log⁡σk−12σk2(yj−μk)2]}\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &=E\left[\log P(y, \gamma \mid \theta) \mid y, \theta^{(i)}\right] \\ &=E\left\{\sum_{k=1}^{K} n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right]\right\} \\ &=\sum_{k=1}^{K}\left\{\sum_{j=1}^{N}\left(E \gamma_{j k}\right) \log \alpha_{k}+\sum_{j=1}^{N}\left(E \gamma_{j k}\right)\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} \end{aligned} Q(θ,θ(i))=E[logP(y,γ∣θ)∣y,θ(i)]=E{k=1∑Knklogαk+j=1∑Nγjk[log(2π1)−logσk−2σk2(yj−μk)2]}=k=1∑K{j=1∑N(Eγjk)logαk+j=1∑N(Eγjk)[log(2π1)−logσk−2σk21(yj−μk)2]}
这里出现了EγjkE \gamma_{j k}Eγjk，依据我们前面的推导它是已知观测数据和当前参数的情况下，隐函数的似然。可以写成：E(γjk∣y,θ)E\left(\gamma_{j k} \mid y, \theta\right)E(γjk∣y,θ)，记为γ^jk\hat{\gamma}_{j k}γ^jk。有
γ^jk=E(γjk∣y,θ)=P(γjk=1∣y,θ)=P(γjk=1,yj∣θ)∑k=1KP(γjk=1,yj∣θ)=P(yj∣γjk=1,θ)P(γjk=1∣θ)∑k=1KP(yj∣γjk=1,θ)P(γjk=1∣θ)=αkϕ(yj∣θk)∑k=1Kαkϕ(yj∣θk),j=1,2,⋯,N;k=1,2,⋯,K\begin{aligned} \hat{\gamma}_{j k} &=E\left(\gamma_{j k} \mid y, \theta\right)=P\left(\gamma_{j k}=1 \mid y, \theta\right) \\ &=\frac{P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)}{\sum_{k=1}^{K} P\left(\gamma_{j k}=1, y_{j} \mid \theta\right)} \\ &=\frac{P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)}{\sum_{k=1}^{K} P\left(y_{j} \mid \gamma_{j k}=1, \theta\right) P\left(\gamma_{j k}=1 \mid \theta\right)} \\ &=\frac{\alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} \mid \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K \end{aligned} γ^jk=E(γjk∣y,θ)=P(γjk=1∣y,θ)=∑k=1KP(γjk=1,yj∣θ)P(γjk=1,yj∣θ)=∑k=1KP(yj∣γjk=1,θ)P(γjk=1∣θ)P(yj∣γjk=1,θ)P(γjk=1∣θ)=∑k=1Kαkϕ(yj∣θk)αkϕ(yj∣θk),j=1,2,⋯,N;k=1,2,⋯,K
推到这可以知道γ^jk\hat{\gamma}_{j k}γ^jk等于当前模型参数下第jjj个观测数据来自第kkk个分模型的概率，称为分模型kkk对观测数据yjy_jyj的响应度。代入QQQ函数可以得到:
Q(θ,θ(i))=∑k=1K{nklog⁡αk+∑j=1Nγ^jk[log⁡(12π)−log⁡σk−12σk2(yj−μk)2]}Q\left(\theta, \theta^{(i)}\right) =\sum_{k=1}^{K}\left\{n_{k} \log \alpha_{k}+\sum_{j=1}^{N}\hat{\gamma}_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} Q(θ,θ(i))=k=1∑K{nklogαk+j=1∑Nγ^jk[log(2π1)−logσk−2σk21(yj−μk)2]}
到此就得到了只含有模型参数的QQQ函数，真的不容易啊！要是没有李航老师的书做参考，估计推到吐血都推不出来！

有了QQQ函数相当于EEE步就有了，MMM步很简单，就是QQQ函数取相应模型参数的偏导然后求其极值点也就是等于0的点即可。求得结果如下：
μ^k=∑j=1Nγ^jkyj∑j=1Nγ^jk,k=1,2,⋯,K\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K μ^k=∑j=1Nγ^jk∑j=1Nγ^jkyj,k=1,2,⋯,K

σ^k2=∑j=1Nγ^jk(yj−μk)2∑j=1Nγ^jk,k=1,2,⋯,K\hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K σ^k2=∑j=1Nγ^jk∑j=1Nγ^jk(yj−μk)2,k=1,2,⋯,K

α^k=nkN=∑j=1Nγ^jkN,k=1,2,⋯,K\hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K α^k=Nnk=N∑j=1Nγ^jk,k=1,2,⋯,K

可以看到其只包含有γ^jk\hat{\gamma}_{j k}γ^jk这一与当前参数相关的变量，因此在实际计算过程中，我们只需要设定初值，然后在EEE步计算出γ^jk\hat{\gamma}_{j k}γ^jk，在M步将计算得到的γ^jk\hat{\gamma}_{j k}γ^jk代入求得新的参数，一直重复直到收敛即可。

真的是推导要老命，编程3分钟啊！

高斯混合模型使用EM算法解决实际问题

已知观测数据 -67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75 试估计两个分量的高斯混合模型的5个参数。

由上面的推导已经知道了所有需要的信息，因此只要设定初值然后直接代公式即可，代码如下：

#include <iostream>
#include <cmath>#define N 15
#define pi 3.1415926535898class theta
{public:double alpha;double mu;double sigma;void print(){std::cout << "--------------" << std::endl;std::cout << "alpha:" << alpha << std::endl;std::cout << "mu:" << mu << std::endl;std::cout << "sigma:" << sigma << std::endl;std::cout << "sigma平方" << sigma* sigma << std::endl;}
};theta mtheta[2];//两个高斯分模型参数
double gamma[2][N];//E步结果gamma
double y[N] = { -67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75 };//观测结果double phi(theta& mtheta, double yj)
{return 1 / (sqrt(2 * pi) * mtheta.sigma) * exp(-pow((yj - mtheta.mu), 2) / (2 * pow(mtheta.sigma, 2)));
}void EStep()
{for (size_t j = 0; j < N; j++){gamma[0][j] = mtheta[0].alpha * phi(mtheta[0], y[j]);gamma[1][j] = mtheta[1].alpha * phi(mtheta[1], y[j]);double sum = gamma[0][j] + gamma[1][j];gamma[0][j] = gamma[0][j] / sum;gamma[1][j] = gamma[1][j] / sum;//std::cout << "gamma0:" << gamma[0][j] << std::endl;//std::cout << "gamma1:" << gamma[1][j] << std::endl;}
}void MStep()
{for (size_t k = 0; k < 2; k++){double mu = 0;double sigma = 0;double sumgamma = 0;for (size_t j = 0; j < N; j++){mu += gamma[k][j] * y[j];sigma += gamma[k][j] * pow((y[j] - mtheta[k].mu), 2);sumgamma += gamma[k][j];}mtheta[k].mu = mu / sumgamma;mtheta[k].sigma = sqrt(sigma / sumgamma);mtheta[k].alpha = sumgamma / N;}
}int main()
{//初始化高斯分模型参数变量mtheta[0].alpha = 0.5;mtheta[0].mu = 30;mtheta[0].sigma = sqrt(500);mtheta[1].alpha = 0.5;mtheta[1].mu = -30;mtheta[1].sigma = sqrt(500);for (size_t k = 0; k < 2; k++){mtheta[k].print();}//迭代10次for (size_t i = 0; i < 10; i++){EStep();MStep();for (size_t k = 0; k < 2; k++){mtheta[k].print();}}
}

如果求出了模型参数后想知道观察结果在这一参数下的分类，再求一次γ^jk\hat{\gamma}_{j k}γ^jk即可！

参考文章

https://zhuanlan.zhihu.com/p/32508410