高斯模型EM算法推导

上文可能花了较多的时间在编辑公式，不过对于它们的由来却没有去深究。

第一眼看到一些结论仍然是很诧异的。虽然和对高斯分布的最大似然估计看起来有几分神似，但是总是觉得不太踏实。此文就是用来梳理推导过程的，在参考文献1中其实有一个手写版的比较详尽的推导，不过仅仅推了μ\mu的计算，而且还有小小的错误。Anyway，先把在Maximization这一步的几个公式再罗列一遍吧

ϕj=1m∑i=1mI{z(i)=j}μj=∑mi=1I{z(i)=j}x(i)∑mi=1I{z(i)=j}Σj=∑mi=1I{z(i)=j}x((i)−μj)(x(i)−μj)T∑mi=1I{z(i)=j}\phi_j=\frac{1}{m}\sum_{i=1}^{m} I\{z^{(i)}=j\} \\ \mu_j=\frac{\sum_{i=1}^{m} I \{z^{(i)}=j\} x^{(i)}}{\sum_{i=1}^{m}I\{z^{(i)}=j\}}\\ \Sigma_j = \frac{\sum_{i=1}^m I\{z^{(i)}=j\}x(^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_{i=1}^{m} I\{z^{(i)}=j\}}

ϕj\phi_j代表着p(z(i)=j)p(z^{(i)}=j)的概率。这个式子是容易理解的，就是所有的m个样本中，输入标签j的数目除以总的样本数就等于ϕj\phi_j.

密度估计的实质

我们进行密度估计的实质，就是在给定m个D维的样本点，x1,...,xm∈RDx_1,...,x_m\in R^D ,令所有的RDR^D上概率密度函数的集合表示成F\mathcal{F},我们要寻找一个概率密度函数f(x)∈Ff(x)\in \mathcal{F},它最有可能产生给定的数据。

一种定义这个函数集F\mathcal{F}的方法就是给所有的函数以相同的数学形式，但是用不同的参数集合θ\theta 来区分它们。比如说，F\mathcal{F}中的函数可以是高斯混合函数：

f(x;θ)=∑k=1KϕkN(x;μk,σk)f(x;\theta)=\sum_{k=1}^{K} \phi_k \mathcal{N}(x; \mu_k,\sigma_k)

由于ff是一个密度函数，所以它必须是非负的，同时积分应该是1。所以我们就可以得到：

1=∫RDf(x;θ)dx=∫RD∑k=1KϕkN(x;μk,σk)dx=11=\int_{R^D} f(x;\theta)dx = \int_{R^D}\sum_{k=1}^{K} \phi_k \mathcal{N}(x;\mu_k,\sigma_k)dx=1

又由于高斯密度函数本身积分为1：

∫RDN(x;μk,Σk)dx=1\int_{R^D} \mathcal{N}(x;\mu_k,\Sigma_k)dx=1

所以可以得到：

∑k=1Kϕk=1\sum_{k=1}^{K} \phi_k=1

所以，我们就能得出这样的结论：

ϕk≥0,∑k=1Kϕk=1\phi_k\ge 0, \sum_{k=1}^{K}\phi_k = 1

这也检验了，我们设定ϕk\phi_k为高斯分量所占的比例是合理的。可以将它称为 mixing probabilities.

推算

使用混合高斯的函数来对cluters进行建模，将每一个cluster分配给一个高斯分量，它的均值靠近cluster的中心，它的标准差就衡量了这个cluster分散的程度。

再来强调一次我们的目的，我们是希望找到这样一个函数，或者说一个模型，通过它最有可能生成我们现在已经获得的数据，所以我们实际上的似然函数最大化就是使得后验概率最大化的过程。

L(x1,...,xm;θ)=∏n=1mf(xn;θ)\mathcal{L}(x_1,...,x_m;\theta)=\prod_{n=1}^{m}f(x_n;\theta)

对于混合高斯函数：

L(x1,...,xm;θ)=∏n=1m∑k=1KϕkN(xn;μk,σk)\mathcal{L(x_1,...,x_m;\theta)}=\prod_{n=1}^{m}\sum_{k=1}^{K}\phi_k \mathcal{N}(x_n; \mu_k,\sigma_k)

对它进行最大化，就是对它的logarithm的最大化，可以表示为：

lnL(x1,...,xm;θ)=∑n=1mln∑k=1KϕkN(xn;μk,σk)ln \mathcal{L}(x1,...,x_m;\theta)=\sum_{n=1}^{m}ln \sum_{k=1}^{K}\phi_k \mathcal{N}(x_n; \mu_k, \sigma_k)

单独对μ\mu求导，为了方便参数仅仅写出μ\mu,可以写成：

ln L(x1,...,xm;μ1,μ2,...μK)=∑n=1mln[∑k=1Kp(xn|wk;μ1,μ2,...μk)p(wk)]ln\ \mathcal{L}(x_1,...,x_m;\mu_1,\mu_2,...\mu_K)=\sum_{n=1}^{m}ln[\sum_{k=1}^{K}p(x_n|w_k;\mu_1,\mu_2,...\mu_k)p(w_k)]

由于：

∂∂xlnf(x)=1f(x)∂f(x)∂x\frac{\partial}{\partial x} ln f(x)=\frac{1}{f(x)}\frac{\partial f(x)}{\partial x}

所以可得：

∂ln L∂μi=∑n=1m1p(xn;μ1,...,μK)∂∂μi∑j=1Kp(xn|wj;μ1,μ2,...,μK)p(wj)=∑n=1m1p(xn;μ1,...,μK)p(wi)∂∂μip(xn|wi;μ1,μ2,...,μK)\frac{\partial ln \ \mathcal{L}}{\partial
\mu_i}=\sum_{n=1}^{m}\frac{1}{p(x_n;\mu_1,...,\mu_K)}\frac{\partial}{\partial \mu_i}\sum_{j=1}^{K}p(x_n|w_j;\mu_1,\mu_2,...,\mu_K)p(w_j) \\=\sum_{n=1}^{m}\frac{1}{p(x_n;\mu_1,...,\mu_K)}p(w_i)\frac{\partial}{\partial \mu_i} p(x_n|w_i;\mu_1,\mu_2,...,\mu_K)

对于高斯正态分布的求导，满足：

∂N(xn;μk,σk)∂μk=N(xn;μk,σk)∂∂μk[−12(||xn−μk||σk)2]\frac{\partial \mathcal{N}(x_n;\mu_k,\sigma_k)}{\partial \mu_k}=\mathcal{N}(x_n;\mu_k,\sigma_k)\frac{\partial}{\partial \mu_k}[-\frac{1}{2}(\frac{||x_n-\mu_k||}{\sigma_k})^2]

并且满足：

∂∂μk||xn−μk||2=∂∂μk(xTnxn+μTkμk−2xTnμk)=2(μk−xn)\frac{\partial}{\partial \mu_k}||x_n-\mu_k||^2=\frac{\partial}{\partial \mu_k}
(x_n^Tx_n+\mu_k^T\mu_k-2x_n^T\mu_k)=2(\mu_k-x_n)

(论文2中的公式写错了吧)

所以上面的式子进一步可以写成：

∂ln L∂μi=∑n=1mp(wi)p(xn|wi;μ1,...,μK)p(xn;μ1,...,μK)(xn−μk)σ2i\frac{\partial ln \ \mathcal{L}}{\partial
\mu_i}=\sum_{n=1}^{m} \frac{p(w_i) p(x_n|w_i;\mu_1,...,\mu_K)}{p(x_n;\mu_1,...,\mu_K)}\frac{(x_n-\mu_k)}{\sigma_i^2}

根据贝叶斯定律，转化成先验概率的形式：

∂ln L∂μi=∑n=1mp(wi|xn;μ1,μ2,...μK)(xn−μk)σ2i\frac{\partial ln \ \mathcal{L}}{\partial
\mu_i}=\sum_{n=1}^{m}p(w_i|x_n;\mu_1,\mu_2,...\mu_K) \frac{(x_n-\mu_k)}{\sigma_i^2}

令偏导数为0，可以得到：

μ^i=∑mn=1P(wi|xn;μ1,...,μK)xn∑mn=1p(wi|xn;μ1,...,μK)\hat\mu_i = \frac{\sum_{n=1}^{m} P(w_i|x_n;\mu_1,...,\mu_K)x_n}{\sum_{n=1}^{m} p(w_i|x_n; \mu_1,...,\mu_K)}

到了这里比较难的就是要继续去推对σ\sigma的微分，为此我试图偷懒从一些文献上找到答案，但是发现大家都很显然的就得到了结果，所以我只能笨笨地自己推了，在参考文献4中貌似有一种更简单的做法。但是当中有个求导的步骤怎么跳过去的我没看懂。

∂Ln L∂σ2=∑n=1mp(wi)p(xn;σ21,...,σ2K)∂∂σ2i[p(xn|wi;σ21,σ22,...,σ2K)]\frac{\partial Ln\ \mathcal{L}}{\partial \sigma^2}=\sum_{n=1}^{m}\frac{p(w_i)}{p(x_n;\sigma_1^2,...,\sigma_K^2)}\frac{\partial}{\partial \sigma_i^2}[p(x_n|w_i;\sigma_1^2,\sigma_2^2,...,\sigma_K^2)]

现在就要回到正态分布的函数形式上：

∂N(xn;μi,σi)∂σi=∂∂σ2[12π−−√σexp(−(xn−μi)22σ2)]=12π−−√σexp(−(xn−μi)22σ2)⋅∂∂σ2(−(xn−μi)22σ2)+exp(−(xn−μi)22σ2)∂∂σ2(12π−−√σ)=N(xn;μi,σi)⋅(xn−μi)22σ4−exp(−(xn−μi)22σ2)⋅122π−−√⋅1σ3=N(xn;μi,σi)⋅(xn−μi)22σ4−N(xn;μi,σi)⋅12σ2\frac{\partial \mathcal{N}(x_n;\mu_i,\sigma_i)}{\partial \sigma_i}=\frac{\partial}{\partial \sigma^2}[\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x_n-\mu_i)^2}{2\sigma^2})]\\=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x_n-\mu_i)^2}{2\sigma^2})\cdot \frac{\partial}{\partial \sigma^2}(-\frac{(x_n-\mu_i)^2}{2\sigma^2})+exp(-\frac{(x_n-\mu_i)^2}{2\sigma^2})\frac{\partial}{\partial \sigma^2}(\frac{1}{\sqrt{2\pi}\sigma})\\=\mathcal{N}(x_n;\mu_i,\sigma_i)\cdot \frac{(x_n-\mu_i)^2}{2\sigma^4}-exp(-\frac{(x_n-\mu_i)^2}{2\sigma^2})\cdot \frac{1}{2\sqrt{2\pi}}\cdot \frac{1}{\sigma^3}\\=
\mathcal{N}(x_n;\mu_i,\sigma_i)\cdot \frac{(x_n-\mu_i)^2}{2\sigma^4}-\mathcal{N}(x_n;\mu_i,\sigma_i)\cdot \frac{1}{2\sigma^2}

这样，上式子变为：

∂Ln L∂σ2i=∑n=1mp(wi)p(xn;σ21,...,σ2K)p(xn|wi,σ21,...,σ2K)[(xn−μi)22σ4i−12σ2i]=∑n=1mp(wi|xn;σ21,σ22,...,σ2K)[(xn−μi)22σ4i−12σ2i]\frac{\partial Ln\ \mathcal{L}}{\partial \sigma_i^2}=\sum_{n=1}^{m}\frac{p(w_i)}{p(x_n;\sigma_1^2,...,\sigma_K^2)}p(x_n|wi,\sigma_1^2,...,\sigma_K^2)[\frac{(x_n-\mu_i)^2}{2\sigma_i^4}-\frac{1}{2\sigma_i^2}]\\
=\sum_{n=1}^{m}p(w_i|x_n;\sigma_1^2,\sigma_2^2,...,\sigma_K^2)[\frac{(x_n-\mu_i)^2}{2\sigma_i^4}-\frac{1}{2\sigma_i^2}]

令偏微分为0，求得结果为：

σ^2i=∑mn=1P(wi|xn;σ21,...,σ2K)(xn−μ^i)2∑mn=1p(wi|xn;σ21,...,σ2K)\hat\sigma_i^2 = \frac{\sum_{n=1}^{m} P(w_i|x_n;\sigma_1^2,...,\sigma_K^2)(x_n-\hat\mu_i)^2}{\sum_{n=1}^{m} p(w_i|x_n; \sigma_1^2,...,\sigma_K^2)}

至于对混合系数的推导，自我感觉这个公式很make sense, 但是实际的推导却要借助拉格朗日乘法因子来引入限制。构造新的目标函数：

J=ln L+λ⋅(1−∑Kj=1wj)J= ln \ \mathcal{L} + \lambda \cdot (1-\sum_{j=1}^{K}w_j)

参考文献:

[1]http://www.cs.cmu.edu/~awm/doc/gmm-algebra.pdf

[2]https://www.cs.duke.edu/courses/spring04/cps196.1/handouts/EM/tomasiEM.pdf

[3]http://www.atmos.washington.edu/~dennis/MatrixCalculus.pdf

[4]http://bengio.abracadoudou.com/lectures/gmm.pdf