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

最新推荐文章于 2024-04-06 23:02:02 发布

_寒潭雁影

最新推荐文章于 2024-04-06 23:02:02 发布

阅读量576

点赞数 2

分类专栏：机器学习数学基础文章标签：高斯混合模型 EM算法

本文链接：https://blog.csdn.net/weixinhum/article/details/108776857

版权

机器学习数学基础专栏收录该内容

11 篇文章 21 订阅

订阅专栏

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

高斯混合模型

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

具体的，高斯混合模型是指具有如下形式的概率分布模型：
$\mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right)$
其中， $\alpha_{k}$ 是系数， $\alpha_{k} \geqslant 0$ , $\quad \sum_{k=1}^{K} \alpha_{k}=1$ ； $\phi\left(y \mid \theta_{k}\right)$ 是高斯分布密度， $\theta_{k}=\left(\mu_{k}, \sigma_{k}^{2}\right)$ ，
$\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)$
称为第 $k$ 个分模型。

$Q$ 函数的一般表达

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

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

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

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

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

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

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

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

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

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

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

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

$\cdots, N ; \quad k=1,2, \cdots, K$

有了观测数据 $y_j$ 及未观测数据 $\gamma_{j k}$ ，那么完全数据是
$\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}\right), \quad j=1,2, \cdots, N$
于是，可以写出完全数据的似然函数：
$\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}$
式中， $n_{k}=\sum_{j=1}^{N} \gamma_{j k}$ ， $\sum_{k=1}^{K} n_{k}=N$ 。

那么，完全数据的对数似然函数为：
$\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]\\$
整个 $Q$ 函数为：
$\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}$
这里出现了 $\gamma_{j k}$ ，依据我们前面的推导它是已知观测数据和当前参数的情况下，隐函数的似然。可以写成： $E\left(\gamma_{j k} \mid y, \theta\right)$ ，记为 $\hat{\gamma}_{j 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}$
推到这可以知道 $\hat{\gamma}_{j k}$ 等于当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $y_j$ 的响应度。代入 $Q$ 函数可以得到:
$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$ 函数，真的不容易啊！要是没有李航老师的书做参考，估计推到吐血都推不出来！

有了 $Q$ 函数相当于 $E$ 步就有了， $M$ 步很简单，就是 $Q$ 函数取相应模型参数的偏导然后求其极值点也就是等于0的点即可。求得结果如下：
$\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$

$\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$

$\hat{\alpha}_{k}=\frac{n_{k}}{N}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K$

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

真的是推导要老命，编程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.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();
		}
	}
}

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

参考文章

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

_寒潭雁影

关注

2
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
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).
复制链接

扫一扫