本篇文章是之前期望极大算法(EM算法)文章的后续，有需要可以先看看那篇文章关于EM算法的推导。
高斯混合模型
高斯混合模型是研究算法的人避不开的一个东西，其在非深度学习的远古时代经常被用到，比如图像处理任务的前背景提取，点云处理任务的点云聚类等等等等。
具体的，高斯混合模型是指具有如下形式的概率分布模型：
称为第k 个分模型。
Q 函数的一般表达
在算法处理的过程中，将问题建模成高斯混合模型后，往往需要去解模型中的参数，这个时候就需要用到EM算法。
在期望极大算法(EM算法)文章的分析中我们已经知道，要进行EM算法，得先得到Q 函数：
在随后的抛硬币问题中我们也介绍了一种求解这个Q QQ函数的方法.但如果看过李老师的书的人会发现，本人没有用书中给出的看起来更具总结性的Q QQ函数形式来求解问题，原因在于当时本人也不明白那个式子是怎么来的。。。
在李航老师的书上，Q 函数是这样定义的：完全数据的对数似然函数log P ( Y , Z ∣ θ )关于在给定观测数据Y YY和当前参数θ ( i )下对未观测数据Z的条件概率分布P ( Z ∣ Y , θ ( i
) ) 的期望称为Q函数，即：
第一次看到这个定义和下面的公式感觉整个人都不好了！实在不知道他们之间为什么是个相等的关系。在经过一两个小时的发呆、无助、掉头发后突然看懂了！
假设log P ( Y , Z ∣ θ ) 是只与变量Z 相关的函数，则可以把其写成f ( Z ) ，当Z取得一个定值的时候，其就是一个固定的数值。如
果这个时候对它取期望，就有：
而如果Z的取值本身受到别的参数x xx，y yy影响，而这些参数都已经给出，则原式可以写成：
代入我们已知的量，等式得证。
由证明的过程也可以知道，这里的给定观测数据Y YY和当前参数两个已知量影响的只有Z ZZ这一变量。由于他们与log P ( Y , Z ∣ θ ) \log P(Y, Z \mid \theta)logP(Y,Z∣θ)中的两个
参数看起来似乎有关系，因此才大大增加了理解的难度。
有了Q QQ函数的一般表达，从式子的形式我们知道关键的一步是把高斯混合模型的log P ( Y , Z ∣ θ ) \log P(Y, Z \mid \theta)logP(Y,Z∣θ)给列出来，也就是把其完全数据的似然函
数的对数列出来。
高斯混合模型参数估计的EM算法
假设数据y 1 , y 2 , ⋯ , y N 由高斯混合模型生成，
其中，θ = ( α 1 , α 2 , ⋯ , α K ; θ 1 , θ 2 , ⋯ , θ K )，求高斯混合模型的参数，就是用EM算法估计高斯混合模型的参数θ 。回顾观测数据y j , j = 1 , 2 , ⋯ , N , 的产生过程：
依概率α k 选择第k 个高斯分布分模型ϕ ( y ∣ θ k ) ；
依第k kk个分模型的概率分布ϕ ( y ∣ θ k ) 生成观测数据y j 。
在高斯混合模型建模中，结果通常是已知的—观测数据y j , j = 1 , 2 , ⋯ , N , 是已知的；
而结果产生自哪个分模型未知—反映观测数据y j 来自第k kk个分模型的数据未知，k = 1 , 2 , ⋯ , K ， k=1,2, 表示，可定义如下:
这里出现了E γ j k ，依据我们前面的推导它是已知观测数据和当前参数的情况下，隐函数的似然。可以写成：E ( γ j k ∣ y , θ ) ，记为γ ^ j k。有
推到这可以知道γ ^ j k 等于当前模型参数下第j 个观测数据来自第k 个分模型的概率，称为分模型k 对观测数据y j 的响应度。代入Q函数可以得到:
到此就得到了只含有模型参数的Q 函数，真的不容易啊！要是没有李航老师的书做参考，估计推到吐血都推不出来！
有了Q 函数相当于E 步就有了，M步很简单，就是Q 函数取相应模型参数的偏导然后求其极值点也就是等于0的点即可。求得结果如下：
可以看到其只包含有γ ^ j k 这一与当前参数相关的变量，因此在实际计算过程中，我们只需要设定初值，然后在E 步计算出γ ^ j k，在M步将计算得到的γ ^ j k 代入求得新的参数，一直重复直到收敛即可。
真的是推导要老命，编程3分钟啊！
高斯混合模型使用EM算法解决实际问题
已知观测数据 -67，-48，6，8，14，16，23，24，28，29，41，49，56，60，75 试估计两个分量的高斯混合模型的5个参数。
由上面的推导已经知道了所有需要的信息，因此只要设定初值然后直接代公式即可，代码如下：
#include
#include
#define N 15
#define pi 3.1415926535898
class 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();
}
}
}如果求出了模型参数后想知道观察结果在这一参数下的分类，再求一次γ ^ j k即可！
参考文章
https://zhuanlan.zhihu.com/p/32508410
