【未完成】混合高斯模型

单一高斯分布情况

在一维条件下使用MLE对未知参数$\mu$和${\sigma }^2$ 进行估计，首先写出
$$p(x)=\prod _{n=1}^{N}N(x_n|\mu ,\sigma )=\frac{1}{(2\pi {\sigma }^2{)}^{\frac{N}{2}}}exp[-\frac{1}{2{\sigma }^2}\sum _{n=1}^N(x_n-\mu {)}^2],$$
取ln后得到：
$$lnp(x)=-\frac{N}{2}ln(2\pi )-\frac{N}{2}ln({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{n=1}^N(x_n-\mu {)}^2,$$
分别对$\mu$和${\sigma }^2$求偏导并令导数为0即可得到估计值。

例如对$\mu$求导得到：
$$\begin{gathered} \frac{\partial lnp(x)}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{n=1}^N(x_n-\mu )=0 \\ \hat{\mu }=\frac{1}{N}\sum _{n=1}^N{x}_n \end{gathered}$$
只与观测值有关。

混合高斯情况

在一维情况下假设存在两个高斯分布，假设方差相同，参数分别是${\mu }_1,{\mu }_2,\sigma$。假设概率分别为p和1-p，可以写出$p(x)$为：
$$p(x)=pN(x|{\mu }_1,\sigma )+(1-p)N(x|{\mu }_2,\sigma )$$
对其取对数后并对${\mu }_1$求偏导得到：
$$\frac{\partial lnp(x)}{\partial {\mu }_1}=\sum _{n=1}^N\frac{pN(x_n|{\mu }_1,\sigma )\frac{x_n-{\mu }_1}{{\sigma }^2}}{pN(x|{\mu }_1,\sigma )+(1-p)N(x|{\mu }_2,\sigma )}$$
令其为0得到：
$${\hat{\mu }}_1=\frac{{\sum }_{n=1}^N\frac{pN(x_n|{\mu }_1,\sigma )x_n}{pN(x|{\mu }_1,\sigma )+(1-p)N(x|{\mu }_2,\sigma )}}{{\sum }_{n=1}^N\frac{pN(x_n|{\mu }_1,\sigma )}{pN(x|{\mu }_1,\sigma )+(1-p)N(x|{\mu }_2,\sigma )}}$$
可见估计值是与带估计值有关的，这使得MLE不能得到解析解。

EM方法来源

如果在上式估计中我们已知${\sum }_{n=1}^N\frac{pN(x_n|{\mu }_1,\sigma )}{pN(x|{\mu }_1,\sigma )+(1-p)N(x|{\mu }_2,\sigma )}$，由贝叶斯定理可以知道这其实就是在观测到$x$时候对应其属于第1类的后验概率$p(z_1|x)$。

更通用的混合模型表达形式为$lnp(x|\theta )$，其中$\theta$为待估计参数，利用全概率公式可以写为：
$$\begin{array}{rl}lnp(x|\theta )& =ln\sum _{z}p(x,z|\theta )\\ & =ln\sum _{z}q(z)\frac{p(x,z|\theta )}{q(z)}\\ & =ln{E}_{z\sim q}[\frac{p(x,z|\theta )}{q(z)}]\end{array}$$
由于上式中包含了log-sum形式不利于求导，利用Jensen不等式可以得到：
$$ln{E}_{z\sim q}[\frac{p(x,z|\theta )}{q(z)}]\ge E_{z\sim q}[ln\frac{p(x,z|\theta )}{q(z)}]=\sum _{z}q(z)ln\frac{p(x,z|\theta )}{q(z)}$$
可以将log-sum形式去除掉，同时最大化$E_{z\sim q}[ln\frac{p(x,z|\theta )}{q(z)}]$可以保证$ln{E}_{z\sim q}[\frac{p(x,z|\theta )}{q(z)}]$也在增大，具体形式如下图所示：

由Jensen不等式知在等号成立时候，$X=E[X]$，此时可以得到：
$$q(z)=p(z|x,\theta )$$
即后验概率。

用EM方法求解GMM

假设存在K种高斯分布混合在一起得到一组观测值$x$，每类高斯分布的均值和方差分别为${\mu }_k$和${\Sigma }_k$
$$p(x)=\prod _{n=1}^{N}p(x_n|\theta )，$$
其中$\theta$表示观测值$x_n$的均值和方差，但是这样任然不能得到正确的分布形式，不知道每一个观测值$x_n$对应的均值和方差，此时数据为incomplete-data。因此需要引入一个参数$z$来表示其属性，得到complete-data。

在概率中，我们观测到的x都是有一个属性的，也就是它的label，完整形式应该是$p(x_n,z_k)$。于是可以改写为：
$$\begin{array}{rl}p(x)& =\sum _{j=1}^{K}p(x,z_j)\\ & =\prod _{n=1}^N\sum _{j=1}^{K}p(x_n,z_j)\\ & =\prod _{n=1}^N\sum _{j=1}^{K}p(z_j)p(x_n|z_j)\\ & =\prod _{n=1}^N\sum _{j=1}^{K}p(z_j)N(x_n|{\mu }_j,{\Sigma }_j)，\end{array}$$

其中$p(z_k)$就表示第k类高斯分布出现的概率，而$p(x_n|z_k)=N(x_n|{\mu }_k,{\Sigma }_k)$，我们就可以写出具体的表达式了：
$$N(x_n|{\mu }_k,{\Sigma }_k)=\frac{1}{\sqrt{(2\pi {)}^{D}det(\Sigma )}}exp[-\frac{1}{2}(x_n-{\mu }_k{)}^T{\Sigma }^{-1}(x_n-{\mu }_k)].$$
再对其取ln后得到：
$$lnp(x)=\sum _{n=1}^{N}ln\sum _{j=1}^K{\pi }_{k}N(x_n|{\mu }_k,{\Sigma }_k)$$

接下来我们不需要使用Jensen不等式，Jensen不等式在GMM问题中仅仅提供了隐函数的概念，并给出具体的求法。

EM（Expectation Maximum）最大期望法，包含E步和M步，分别是求出后验概率（期望），然后通过期望来优化参数。其中E步和MLE一样需要对似然函数求偏导并令其为0，需要用到很多矩阵求导方法，可以在这个网站https://www.matrixcalculus.org/上得到结果。

E步

对${\mu }_k$求偏导并令其为0得到：
$$\frac{\partial lnp(x)}{\partial {\mu }_k}=\sum _{n=1}^N\frac{{\pi }_{k}N(x_n|{\mu }_k,{\Sigma }_k)(-{\Sigma }_k^{-1})({\mu }_k-x_n)}{{\sum }_{j=1}^K{\pi }_{j}N(x_n|{\mu }_j,{\Sigma }_j)}$$
$$\sum _{n=1}^N\gamma (z_{nk})({\mu }_k-x_n)=0$$
$${\hat{\mu }}_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})x_n$$
其中$p(z_k|x_n)=\gamma (z_{nk})=\frac{{\pi }_{k}N(x_n|{\mu }_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_{k}N(x_n|{\mu }_j,{\Sigma }_j)}$表示观测到$x_n$属于为$z_k$类的后验概率，也可以叫作responsibility，$N_k={\sum }_{n=1}^N\gamma (z_{nk})$表示k类的个数。因此每一个观测到的$x_n$乘上权重其属于$z_k$类的概率权重$\gamma (z_{nk})$之后再除以k类的总数$N_k$，就得到了第k类的均值${\mu }_k$。

对${\Sigma }_k$求导并令其为0得到：
$$\frac{\partial lnp(x)}{\partial {\Sigma }_k}=\sum _{n=1}^N\gamma (z_{nk})(-\frac{1}{2})({\Sigma }^{-1}{\mu }_k-x_n)({\mu }_k-x_n{)}^T{\Sigma }^{-1})$$
$${\hat{\Sigma }}_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})({\mu }_k-x_n)({\mu }_k-x_n{)}^T.$$
该公式也是符合协方差定义$\Sigma =E[x{x}^T]$，并在此基础上进行加权从而得到k类的协方差。

求$\pi$时因为有一个约束，因此我们需要用拉格朗日法引入约束得到目标函数:
$$lnp(x)+\lambda (\sum _{j=1}^K{\pi }_j-1),$$
在此基础上对${\pi }_k$求导，并令其为0得到：
$$\frac{\partial [lnp(x)+\lambda ({\sum }_{j=1}^K{\pi }_j-1)]}{\partial {\Sigma }_k}=\sum _{n=1}^N\frac{N(x_n|{\mu }_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_{j}N(x_n|{\mu }_j,{\Sigma }_j)}+\lambda$$
两边同时乘上${\pi }_j$并对$\pi$进行累加或者将$\lambda$代入约束可以得到：
$$\sum _{j=1}^K{\pi }_j\sum _{n=1}^N\frac{N(x_n|{\mu }_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_{j}N(x_n|{\mu }_j,{\Sigma }_j)}+\sum _{j=1}^K\lambda {\pi }_j=0$$
$$\lambda =-N$$
代入$\lambda$后得到：
$${\pi }_k=\frac{N_k}{N}，$$
显然这个公式的物理意义是和先前定义${\pi }_k$的物理意义一致的。

可以看出每一个参数都与$\gamma (z_{nk})=p(z_k|x_n)$相关的，因此只要在E步求出最新的这个值，就可以在M步中更新三个参数了。

M步

在E布我们得到了以下三个参数的迭代方式：
$${\hat{\mu }}_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})x_n$$
$${\hat{\Sigma }}_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})({\mu }_k-x_n)({\mu }_k-x_n{)}^T$$
$${\pi }_k=\frac{N_k}{N}$$
后就可以根据迭代法来更新参数，最终收敛时得到估计值。

以下是来自参考书的步骤：

Matlab代码

待完善

参考资料

《Pattern Recognition and Machine Learning》—— Christopher M. Bishop
【机器学习】EM——期望最大（非常详细）
Gaussian mixture models and the EM algorithm —— Ramesh Sridharan