GMM高斯混合模型学习笔记(EM算法求解)
2015-06-11 23:02
337 查看
提出混合模型主要是为了能更好地近似一些较复杂的样本分布,通过不断增加component个数,可以任意地逼近任何连续的概率分布,所以我们认为任何样本分布都可以用混合模型来建模。因为高斯函数具有一些很实用的性质,所以高斯混合模型被广泛地使用。
GMM与kmeans类似,也是属于clustering,不同的是,kmeans是把每个样本点聚到其中一个cluster,而GMM是给出这些样本点到每个cluster的概率,每个component就是一个聚类中心。
GMM(Gaussian Mixture Model)高斯混合模型,由K个不同的Gaussian线性组合而成,每个Gaussian是混合模型的一个component,GMM的概率密度函数如下:
p(x)=∑k=1Kp(k)(x|k)=∑k=1Kπk(x|μk,∑k)p(x)=\sum_{k=1}^Kp(k)(x|k) =\sum_{k=1}^K\pi_k\mathcal N(x|\mu_k,\sum_k)
根据上式,从GMM中生成一个样本点x分两步:
1,从K个component中随机的选择一个
2,从该component中选择一个点
参数说明:N个样本点,K个component,μk,∑k\mu_k,\sum_k 是第k个component的均值和协方差矩阵,是模型参数,是需要估计的。πk\pi_k是mixing coefficient,表示第k个component被选中的概率,πk=1N∑Nn=1znk\pi_k=\frac{1}{N}\sum_{n=1}^N\mathbf z_{nk},也是模型参数,需要估计。N是高斯(正态)分布。
对一个样本集建立高斯混合模型的过程,就是根据已知样本集X反推高斯混合模型的参数(μ,∑,π\mu,\sum,\pi),这是一个参数估计问题。首先想到用最大似然的方法求解,也就是,要确定参数π,μ,∑\pi,\mu,\sum使得它所确定的概率分布生成这些样本点的概率最大,这个概率也就是似然函数,如下:
p(x)=∏n=1Np(xi)p(x)=\prod_{n=1}^Np(x_i)
而一般对于单个样本点其概率较小,多个相乘后更小,容易造成浮点数下溢,所以一般是对似然函数求log,变成加和形式:
∑i=1Nlnp(xi)\sum_{i=1}^Nlnp(x_i)
这个叫做log似然函数,目标是要最大化它。用log似然函数对参数分别求偏导,令偏导等于0,可求解得参数。
然而,GMM的log似然函数是如下形式:
lnp(X)=∑i=1Nln[∑k=1Kπk(xi|μk,∑k)]lnp(X)=\sum_{i=1}^Nln[\sum_{k=1}^K\pi_k\mathcal N(x_i|\mu_k,\sum_k)]
可以看到对数中有求和,直接求导求解将导致一系列复杂的运算,故考虑使用EM算法。(具体思路见上一篇:EM算法学习笔记)
考虑GMM生成一个样本点的过程,这里对每个xi\mathbf x_i引入隐变量z,z是一个K维向量,如果生成xi\mathbf x_i时选择了第k个component,则zk=1\mathbf z_k=1,其他元素都为0,∑Kk=1zk=1\sum_{k=1}^K\mathbf z_k=1.
假设z是已知的,则样本集变成了{X,Z},要求解的似然函数变成了:
p(X,Z|μ,∑,π)=∏n=1N∏k=1Kπznkk(xn|μk,∑k)znkp(X,Z|\mu,\sum,\pi)=\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}\mathcal N(\mathbf x_n|\mu_k,\sum_k)^{z_{nk}}
log似然函数为:
lnp(X,Z|μ,∑,π)=∑n=1N∑k=1Kznk[lnπk+ln(xn|μk,∑k)].(∗)lnp(X,Z|\mu,\sum,\pi)=\sum_{n=1}^N\sum_{k=1}^K\mathbf z_{nk}[ln\pi_k + ln\mathcal N(\mathbf x_n|\mu_k,\sum_k)].(*)
可以看到,这次ln直接对Gaussian作用,求和在ln外面,所以可以直接求最大似然解了。
1,初始化一组模型参数π,μ,∑\pi,\mu,\sum
2,E-step
然而,事实上z是不知道的,我们只是假设z已知。而z的值是通过后验概率观测,所以这里考虑用z值的期望在上述似然函数中代替z。
对于一个样本点x\mathbf x:
p(z)=∏k=1Kπzkkp(\mathbf z)=\prod_{k=1}^K\pi_k^{z_k}
p(x|zk=1)=(x|μk,∑k)p(\mathbf x|\mathbf z_k=1)=\mathcal N(x|\mu_k,\sum_k)
p(x|z)=∏k=1K(x|μk,∑k)zkp(\mathbf x|\mathbf z)=\prod_{k=1}^K\mathcal N(\mathbf x|\mu_k,\sum_k)^{z_k}
p(x)=∑zp(z)p(x|z)=∑k=1Kπk(x|μk,∑k)p(\mathbf x)=\sum_zp(\mathbf z)p(\mathbf x|\mathbf z)=\sum_{k=1}^K\pi_k\mathcal N(\mathbf x|\mu_k,\sum_k)
后验概率(固定μ,∑,π\mu,\sum,\pi):
p(z|x,μ,∑,π)=p(x|z)p(z)p(x)正比于∏n=1N∏k=1K[πk(xn|μk,∑k)]znkp(\mathbf z|\mathbf x,\mu,\sum,\pi)=\frac{p(\mathbf x|\mathbf z)p(\mathbf z)}{p(\mathbf x)}正比于\prod_{n=1}^N\prod_{k=1}^K[{\pi_k\mathcal N(x_n|\mu_k,\sum_k)}]^{z_{nk}}
因为{zn\mathbf z_n}之间是相互独立的。
计算z期望γ(znk)\gamma(\mathbf z_{nk})(z向量只有一个值取1,其余为0):
γ(znk)=E[znk]=0∗p(znk=0|xn)+1∗p(znk=1|xn)=p(znk=1|xn)=p(znk=1)p(xn|znk=1)p(xn)=πk(x|μk,∑k)∑Kj=1πj(x|μj,∑j).\gamma(\mathbf z_{nk})=E[\mathbf z_{nk}]=0*p(\mathbf z_{nk}=0|\mathbf x_n)+1*p(\mathbf z_{nk}=1|\mathbf x_n)=p(\mathbf z_{nk}=1|\mathbf x_n)=\frac{p(\mathbf z_{nk}=1)p(\mathbf x_n|\mathbf z_{nk}=1)}{p(\mathbf x_n)}=\frac{\pi_k\mathcal N(\mathbf x|\mu_k,\sum_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x|\mu_j,\sum_j)}.
将z值用期望代替,则待求解的log似然函数(*)式变为:
Ez[lnp(X,Z|μ,∑,π)]=∑n=1N∑k=1Kγ(znk)[lnπk+ln(xn|μk,∑k)].E_z[lnp(X,Z|\mu,\sum,\pi)]=\sum_{n=1}^N\sum_{k=1}^K\gamma (\mathbf z_{nk})[ln\pi_k + ln\mathcal N(\mathbf x_n|\mu_k,\sum_k)].
3,M-step
现在可以最大化似然函数求解参数了,首先对μ\mu求偏导,令偏导等于0,可得:
∑n=1N∑k=1Kγ(znk)∑k(xn−μk)=0\sum_{n=1}^N\sum_{k=1}^K\gamma (\mathbf z_{nk})\sum_k(\mathbf x_n-\mu_k)=0
μk=1Nk∑n=1Nγ(znk)xn,其中Nk=∑n=1Nγ(znk).\mu_k=\frac{1}{N_k}\sum_{n=1}^N\gamma (\mathbf z_{nk}){\mathbf x_n},其中N_k=\sum_{n=1}^N\gamma (\mathbf z_{nk}).
NkN_k 是“the effective number of points assigned to cluster k”.
再对∑k\sum_k求偏导,令偏导等于0,可得:
∑k=1Nk∑n=1Nγ(znk)(xn−μk)(xn−μk)T\sum_k=\frac{1}{N_k}\sum_{n=1}^N\gamma (\mathbf z_{nk})(\mathbf x_n-\mu_k)(\mathbf x_n-\mu_k)^T
接下来还需求解π\pi,注意到π\pi需满足∑Kk=1πk=1\sum_{k=1}^K\pi_k=1,所以这是一个带等式约束的最大值问题,使用拉格朗日乘数法。
构造拉格朗日函数:
L=lnp(X|π,μ,∑)+λ(∑k=1Kπk−1).L=lnp(X|\pi,\mu,\sum)+\lambda(\sum_{k=1}^K\pi_k-1).
对π\pi求导,令导数为0:
∑n=1N(x|μk,∑k)∑Kj=1πj(x|μj,∑j)+λ=0\sum_{n=1}^N\frac{\mathcal N(\mathbf x|\mu_k,\sum_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x|\mu_j,\sum_j)}+\lambda=0
两边同乘πk\pi_k得:
∑n=1Nγ(znk)+λπk=0\sum_{n=1}^N\gamma (\mathbf z_{nk}) + \lambda\pi_k=0
Nk+λπk=0N_k+\lambda\pi_k=0
两边对k求和:
∑k=1KNk+∑k=1Kλπk=0\sum_{k=1}^KN_k+\sum_{k=1}^K\lambda\pi_k=0
N+λ=0N+\lambda=0
可得:λ=−N\lambda=-N
代入可得:πk=NkN.\pi_k=\frac{N_k}{N}.
4,检查是否收敛
重复E-step和M-step两步,直到收敛,即可求得一个局部最优解。
GMM的建模过程如下图(k=2,高斯分布是蓝色和红色圈):
主要参考资料:
《Pattern Recognization and Machine Learning》
帮助理解:
http://blog.pluskid.org/?p=39
GMM与kmeans类似,也是属于clustering,不同的是,kmeans是把每个样本点聚到其中一个cluster,而GMM是给出这些样本点到每个cluster的概率,每个component就是一个聚类中心。
GMM(Gaussian Mixture Model)高斯混合模型,由K个不同的Gaussian线性组合而成,每个Gaussian是混合模型的一个component,GMM的概率密度函数如下:
p(x)=∑k=1Kp(k)(x|k)=∑k=1Kπk(x|μk,∑k)p(x)=\sum_{k=1}^Kp(k)(x|k) =\sum_{k=1}^K\pi_k\mathcal N(x|\mu_k,\sum_k)
根据上式,从GMM中生成一个样本点x分两步:
1,从K个component中随机的选择一个
2,从该component中选择一个点
参数说明:N个样本点,K个component,μk,∑k\mu_k,\sum_k 是第k个component的均值和协方差矩阵,是模型参数,是需要估计的。πk\pi_k是mixing coefficient,表示第k个component被选中的概率,πk=1N∑Nn=1znk\pi_k=\frac{1}{N}\sum_{n=1}^N\mathbf z_{nk},也是模型参数,需要估计。N是高斯(正态)分布。
对一个样本集建立高斯混合模型的过程,就是根据已知样本集X反推高斯混合模型的参数(μ,∑,π\mu,\sum,\pi),这是一个参数估计问题。首先想到用最大似然的方法求解,也就是,要确定参数π,μ,∑\pi,\mu,\sum使得它所确定的概率分布生成这些样本点的概率最大,这个概率也就是似然函数,如下:
p(x)=∏n=1Np(xi)p(x)=\prod_{n=1}^Np(x_i)
而一般对于单个样本点其概率较小,多个相乘后更小,容易造成浮点数下溢,所以一般是对似然函数求log,变成加和形式:
∑i=1Nlnp(xi)\sum_{i=1}^Nlnp(x_i)
这个叫做log似然函数,目标是要最大化它。用log似然函数对参数分别求偏导,令偏导等于0,可求解得参数。
然而,GMM的log似然函数是如下形式:
lnp(X)=∑i=1Nln[∑k=1Kπk(xi|μk,∑k)]lnp(X)=\sum_{i=1}^Nln[\sum_{k=1}^K\pi_k\mathcal N(x_i|\mu_k,\sum_k)]
可以看到对数中有求和,直接求导求解将导致一系列复杂的运算,故考虑使用EM算法。(具体思路见上一篇:EM算法学习笔记)
考虑GMM生成一个样本点的过程,这里对每个xi\mathbf x_i引入隐变量z,z是一个K维向量,如果生成xi\mathbf x_i时选择了第k个component,则zk=1\mathbf z_k=1,其他元素都为0,∑Kk=1zk=1\sum_{k=1}^K\mathbf z_k=1.
假设z是已知的,则样本集变成了{X,Z},要求解的似然函数变成了:
p(X,Z|μ,∑,π)=∏n=1N∏k=1Kπznkk(xn|μk,∑k)znkp(X,Z|\mu,\sum,\pi)=\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}\mathcal N(\mathbf x_n|\mu_k,\sum_k)^{z_{nk}}
log似然函数为:
lnp(X,Z|μ,∑,π)=∑n=1N∑k=1Kznk[lnπk+ln(xn|μk,∑k)].(∗)lnp(X,Z|\mu,\sum,\pi)=\sum_{n=1}^N\sum_{k=1}^K\mathbf z_{nk}[ln\pi_k + ln\mathcal N(\mathbf x_n|\mu_k,\sum_k)].(*)
可以看到,这次ln直接对Gaussian作用,求和在ln外面,所以可以直接求最大似然解了。
1,初始化一组模型参数π,μ,∑\pi,\mu,\sum
2,E-step
然而,事实上z是不知道的,我们只是假设z已知。而z的值是通过后验概率观测,所以这里考虑用z值的期望在上述似然函数中代替z。
对于一个样本点x\mathbf x:
p(z)=∏k=1Kπzkkp(\mathbf z)=\prod_{k=1}^K\pi_k^{z_k}
p(x|zk=1)=(x|μk,∑k)p(\mathbf x|\mathbf z_k=1)=\mathcal N(x|\mu_k,\sum_k)
p(x|z)=∏k=1K(x|μk,∑k)zkp(\mathbf x|\mathbf z)=\prod_{k=1}^K\mathcal N(\mathbf x|\mu_k,\sum_k)^{z_k}
p(x)=∑zp(z)p(x|z)=∑k=1Kπk(x|μk,∑k)p(\mathbf x)=\sum_zp(\mathbf z)p(\mathbf x|\mathbf z)=\sum_{k=1}^K\pi_k\mathcal N(\mathbf x|\mu_k,\sum_k)
后验概率(固定μ,∑,π\mu,\sum,\pi):
p(z|x,μ,∑,π)=p(x|z)p(z)p(x)正比于∏n=1N∏k=1K[πk(xn|μk,∑k)]znkp(\mathbf z|\mathbf x,\mu,\sum,\pi)=\frac{p(\mathbf x|\mathbf z)p(\mathbf z)}{p(\mathbf x)}正比于\prod_{n=1}^N\prod_{k=1}^K[{\pi_k\mathcal N(x_n|\mu_k,\sum_k)}]^{z_{nk}}
因为{zn\mathbf z_n}之间是相互独立的。
计算z期望γ(znk)\gamma(\mathbf z_{nk})(z向量只有一个值取1,其余为0):
γ(znk)=E[znk]=0∗p(znk=0|xn)+1∗p(znk=1|xn)=p(znk=1|xn)=p(znk=1)p(xn|znk=1)p(xn)=πk(x|μk,∑k)∑Kj=1πj(x|μj,∑j).\gamma(\mathbf z_{nk})=E[\mathbf z_{nk}]=0*p(\mathbf z_{nk}=0|\mathbf x_n)+1*p(\mathbf z_{nk}=1|\mathbf x_n)=p(\mathbf z_{nk}=1|\mathbf x_n)=\frac{p(\mathbf z_{nk}=1)p(\mathbf x_n|\mathbf z_{nk}=1)}{p(\mathbf x_n)}=\frac{\pi_k\mathcal N(\mathbf x|\mu_k,\sum_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x|\mu_j,\sum_j)}.
将z值用期望代替,则待求解的log似然函数(*)式变为:
Ez[lnp(X,Z|μ,∑,π)]=∑n=1N∑k=1Kγ(znk)[lnπk+ln(xn|μk,∑k)].E_z[lnp(X,Z|\mu,\sum,\pi)]=\sum_{n=1}^N\sum_{k=1}^K\gamma (\mathbf z_{nk})[ln\pi_k + ln\mathcal N(\mathbf x_n|\mu_k,\sum_k)].
3,M-step
现在可以最大化似然函数求解参数了,首先对μ\mu求偏导,令偏导等于0,可得:
∑n=1N∑k=1Kγ(znk)∑k(xn−μk)=0\sum_{n=1}^N\sum_{k=1}^K\gamma (\mathbf z_{nk})\sum_k(\mathbf x_n-\mu_k)=0
μk=1Nk∑n=1Nγ(znk)xn,其中Nk=∑n=1Nγ(znk).\mu_k=\frac{1}{N_k}\sum_{n=1}^N\gamma (\mathbf z_{nk}){\mathbf x_n},其中N_k=\sum_{n=1}^N\gamma (\mathbf z_{nk}).
NkN_k 是“the effective number of points assigned to cluster k”.
再对∑k\sum_k求偏导,令偏导等于0,可得:
∑k=1Nk∑n=1Nγ(znk)(xn−μk)(xn−μk)T\sum_k=\frac{1}{N_k}\sum_{n=1}^N\gamma (\mathbf z_{nk})(\mathbf x_n-\mu_k)(\mathbf x_n-\mu_k)^T
接下来还需求解π\pi,注意到π\pi需满足∑Kk=1πk=1\sum_{k=1}^K\pi_k=1,所以这是一个带等式约束的最大值问题,使用拉格朗日乘数法。
构造拉格朗日函数:
L=lnp(X|π,μ,∑)+λ(∑k=1Kπk−1).L=lnp(X|\pi,\mu,\sum)+\lambda(\sum_{k=1}^K\pi_k-1).
对π\pi求导,令导数为0:
∑n=1N(x|μk,∑k)∑Kj=1πj(x|μj,∑j)+λ=0\sum_{n=1}^N\frac{\mathcal N(\mathbf x|\mu_k,\sum_k)}{\sum_{j=1}^K\pi_j\mathcal N(\mathbf x|\mu_j,\sum_j)}+\lambda=0
两边同乘πk\pi_k得:
∑n=1Nγ(znk)+λπk=0\sum_{n=1}^N\gamma (\mathbf z_{nk}) + \lambda\pi_k=0
Nk+λπk=0N_k+\lambda\pi_k=0
两边对k求和:
∑k=1KNk+∑k=1Kλπk=0\sum_{k=1}^KN_k+\sum_{k=1}^K\lambda\pi_k=0
N+λ=0N+\lambda=0
可得:λ=−N\lambda=-N
代入可得:πk=NkN.\pi_k=\frac{N_k}{N}.
4,检查是否收敛
重复E-step和M-step两步,直到收敛,即可求得一个局部最优解。
GMM的建模过程如下图(k=2,高斯分布是蓝色和红色圈):
主要参考资料:
《Pattern Recognization and Machine Learning》
帮助理解:
http://blog.pluskid.org/?p=39
相关文章推荐
- 2015061109 - 如何管理自己的人脉?
- cocos2d-js添加百度appx的插屏广告(通过jsb反射机制)
- Arduino 入门程序示例之一个 LED(2015-06-11)
- sqlite使用之模糊查询数据库数据的三种方式(待完善)
- [leetcode] Reorder List
- response.setContentType与 request.setCharacterEncoding 区别
- win7下VC6.0出现Unable to register this add-in because its DLLRegisterServer returnan error
- 第五章循环
- 树莓派远程桌面
- 解决HttpServletResponse输出的中文乱码问题
- [软件使用]在Word中使用通配符查询
- 自己设定Macbook风扇转速,让苹果不再发烫的秘笈
- hdoj 3501 【欧拉函数 求小于或者等于n的数中 与n互质的数总和】
- 简单测试补码、类型最大值等
- 引用和指针的区别
- 2015061108 - 离职前一定要找好下家吗?
- 习题3-41
- Leetcode[101]-Symmetric Tree
- 【转】gcc warning: braces around scalar initializer (标量初始化的括号)
- ubuntu13 安装 pycharm