第十五章 EM期望极大算法及其推广

导读

• 概率模型有时既含有观测变量，又含有隐变量或潜在变量。这句很重要，有时候我们能拿到最后的观测结果，却拿不到中间的过程记录

• EM算法可以用于生成模型的非监督学习，EM算法是个一般方法，不具有具体模型

  EM算法是一种迭代算法，用于含有隐变量的概率模型的极大似然估计，或极大后验概率估计

• 似然和概率的关系可以推广了解，这章关于概率和似然的符号表示，概率和似然是同样的形式描述的都是可能性，$P(Y|\theta )$是一个两变量的函数，似然是给定结果，求参数可能性；概率是给定参数求结果可能性
  是不是可以认为，似然和概率分别对应了Train和Predict的过程？

  假设给定观测数据Y，其概率分布是P(Y|θ)，其中θ是概率分布的参数。那么不完全数据Y的似然函数是通过将Y视为已知观测数据，而将未观测到的数据表示为缺失或隐变量。似然函数可以表示为P(Y|θ)

• 无论是三硬币还是GMM，采样的过程都是如下:

  1. Sample $z_i\sim p(z|\pi )$
  2. Sample $x_i\sim p(x|\pi )$

  注意，这里用到了$\pi$，在强化学习中，随机性策略$\pi (x,a)$表示为状态$x$下选择动作$a$的概率。

• 关于EM算法的解释
  注意这里EM不是模型，是个一般方法，不具有具体的模型，这点前面已经提到

  1. PRML
     $kmeans\to GMM\to EM$
     k均值聚类算法就是一个典型的EM算法
  2. 统计学习方法
     1. $MLE\to B$
     2. $F$函数的极大-极大算法

• HMM作了两个基本假设，实际上是在说在图模型中，存在哪些边

符号说明

一般地，用$Y$表示观测随机变量的数据，$Z$表示隐随机变量的数据。$Y$和$Z$一起称为完全数据(complete-data)，观测数据$Y$又称为不完全数据(incomplete-data)

上面这个概念很重要，Dempster在1977年提出EM算法的时候文章题目就是《Maximum likelihood from incomplete data via the EM algorithm》，具体看书中本章参考文献[^3]

假设给定观测数据$Y$，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数
那么不完全数据$Y$的似然函数是$P(Y|\theta )$，对数似然函数是$L(\theta )=\log P(Y|\theta )$

假设$Y$和$Z$的联合概率分布是$P(Y,Z|\theta )$，那么完全数据的对数似然函数是$\log P(Y,Z|\theta )$

上面这部分简单对应一下，这里说明一下，你看到下面概率分布和似然函数形式看起来一样。在概率中，$\theta$已知, 求$Y$，在似然函数中通过已知的Y去求$\theta$

有一点要注意下, 这里没有出现$X$, 有提到一种理解

• 有时训练数据只有输入没有对应的输出$(x_1,\cdot ),(x_2,\cdot ),\dots ,(x_N,\cdot )$，从这样的数据学习模型称为非监督学习问题。
• EM算法可以用于生成模型的非监督学习。
• 生成模型由联合概率分布$P(X,Y)$表示，可以认为非监督学习训练数据是联合概率分布产生的数据。$X$为观测数据， $Y$为未观测数据

有时候，只观测显变量看不到关系，就需要把隐变量引进来。

混合模型

书中用三硬币模型做为引子，在学习这部分内容的时候，注意体会观测数据的作用。

伯努利混合模型(三硬币模型)

问题描述

问题的描述过程中有这样一句：独立的重复$n$次实验(这里$n=10$)，观测结果如下:

1,1,0,1,0,0,1,0,1,1

思考上面这个观测和1,1,1,1,1,1,0,0,0,0有区别么?

没有任何信息的前提下，我们得到上面的观测数据可以假定是一个二项分布的形式，参数$n=10,p=0.6$

把$k=6$次成功分布在$n=10$次试验中有$C(10,6)$种可能

所以上面两个观测序列，可能出自同一个模型。在这个问题的求解上是没有区别的

我们通过一段代码来生成这个数据

import numpy as npp = 0.6
n = 10
# np.random.seed(2018)
flag_a = 1
flag_b = 1
cnt = 0
while flag_a or flag_b:tmp = np.random.binomial(1, p, n)if (tmp == np.array([1,1,1,1,1,1,0,0,0,0])).all():flag_a = 0print("[1,1,1,1,1,1,0,0,0,0] at %d\n" % cnt)if (tmp == np.array([1,1,0,1,0,0,1,0,1,1])).all():flag_b = 0print("[1,1,0,1,0,0,1,0,1,1] at %d\n" % cnt)cnt += 1

实际上题目的描述中说明了观测数据生成的过程，这些参数是未知的，所以需要对这些参数进行估计

三硬币模型可以写作
$$\begin{array}{c}\begin{array}{rl}P(y|\theta )& =\sum _{z}P(y,z|\theta )\\ & =\sum _{z}P(z|\theta )P(y|z,\theta )\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}\end{array}$$
以上

1. 随机变量$y$是观测变量，表示一次试验观测的结果是1或0
2. 随机变量$z$是隐变量，表示未观测到的掷硬币$A$的结果
3. $\theta =(\pi ,p,q)$是模型参数
4. 这个模型是以上数据(1,1,0,1,0,0,1,0,1,1)的生成模型

观测数据表示为$Y=(Y_1,Y_2,Y_3,\dots ,Y_n{)}^T$, 未观测数据表示为$Z=(Z_1,Z_2,Z_3,\dots ,Z_n{)}^T$, 则观测数据的似然函数为

其实觉得这里应该是小写的$y=(y_1,y_2,\dots ,y_n),z=(z_1,z_2,\dots ,z_n)$
$$P(Y|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )$$
注意这里的求和是下面的"+"描述的部分

即
$$P(Y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$

$$p(x|\theta )p(\theta )=p(\theta |x)p(x)\Rightarrow p(\theta |x)=\frac{p(x|\theta )p(\theta )}{p(x)}\Rightarrow \mathrm{argmax}p(\theta |x)$$

极大似然估计：
$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )$$
极大后验概率：
$$\text{argmax}p(\theta |x)⟺\text{argmax}p(x|\theta )p(\theta )$$
均匀分布下，后一部分是等概率的，N

KL散度=交叉熵+const

三硬币模型的EM算法

1.初值

EM算法首选参数初值，记作${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$, 然后迭代计算参数的估计值。

如果第$i$次迭代的模型参数估计值为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$

2.E步

那么第$i+1$ 次迭代的模型参数估计值表示为
$${\mu }_j^{i+1}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_j}}$$
因为是硬币，只有0，1两种可能，所有有上面的表达

这个表达方式还可以拆成如下形式
$${\mu }_j^{i+1}=\{\begin{array}{ll}\frac{{\pi }^{(i)}{p}^{(i)}}{{\pi }^{(i)}{p}^{(i)}+(1-{\pi }^{(i)})q^{(i)}}& ,y_j=1\\ \frac{{\pi }^{(i)}(1-p^{(i)})}{{\pi }^{(i)}(1-p^{(i)})+(1-{\pi }^{(i)})(1-q^{(i)})}& ,y_j=0\end{array}$$
所以, 这步(求${\mu }_j$)干了什么，样本起到了什么作用？

这一步，通过假设的参数，计算了不同的样本对假设模型的响应(${\mu }_j$)，注意这里因为样本($y_j$)是二值的，所以，用 { y j , 1 − y j } \{y_j, 1-y_j\} {yj​,1−yj​} 构成了one-hot的编码，用来表示样本归属的假设

以上，有点绕

这一步是什么的期望？书中有写，**观测数据来自硬币$B$的概率, 在二项分布的情况下, 响应度和概率是一个概念. **这个说明，有助于后面M步公式的理解。

3.M步

$$\begin{array}{rl}{\pi }^{(i+1)}& =\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}\\ p^{(i+1)}& =\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}\\ q^{(i+1)}& =\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})}\end{array}$$
上面，红色部分的公式从观测数据是来自硬币B的概率这句来理解

初值影响

这个例子里面0.5是个合理又牛逼的初值。迭代收敛的最后结果是(0.5, 0.6, 0.6)

这个结果说明，如果A是均匀的，那么一个合理的解就是B，C是同质的。他们的分布情况和观测的分布一致

p,q 含义

这里面$p$对应了$A=1$，$B=1$，$q$对应了$A=0$，$C=1$

这三个公式可以改写成如下形式:
$$\begin{array}{rl}{\pi }^{(i+1)}& =\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}\\ p^{(i+1)}& =\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n({\mu }_j^{(i+1)}{y}_j+{\mu }_j^{(i+1)}(1-y_j)}\\ q^{(i+1)}& =\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n((1-{\mu }_j^{(i+1)})y_j+(1-{\mu }_j^{(i+1)})(1-y_j))}\end{array}$$
$\pi$的表达式回答这样一个问题: 刷了这么多样本，拿到一堆数，那么$\pi$应该是多少，均值是个比较好的选择

$p$的表达式回答这样一个问题: 如果我知道每个结果$y_j$以${\mu }_j$的可能来自硬币B(A=1)，那么用这些数据刷出来他可能是正面的概率。这里面${\mu }_j$对应了$A=1$

$q$的表达式同理，其中$1-{\mu }_j$对应了$A=0$

到后面讲高斯混合模型的时候，可以重新审视这里
$$\begin{array}{rl}{\alpha }_0& \leftrightarrow \pi \\ {\mu }_0& \leftrightarrow p^{y_j}(1-p{)}^{1-y_j}\\ {\alpha }_1& \leftrightarrow 1-\pi \\ {\mu }_1& \leftrightarrow q^{y_j}(1-q{)}^{1-y_j}\end{array}$$
以上对应了包含两个分量的伯努利混合模型

bmm.py对伯努利混合模型做了实现, 有几点说明一下:

1. $(p^{(i)}{)}^{y_i}(1-p^{(i)}{)}^{1-y_i}$这个表达式对应了伯努利分布的概率密度，可以表示成矩阵乘法，尽量不要用for，效率会差

2. 采用了$\pi ,p,q$来表示，注意在题目的说明部分有说明三个符号的含义

3. 实际上不怎么抛硬币，但是0-1的伯努利分布很多

   当参数$\theta$的维数为$d(d\ge 2)$的时候，可以采用一种特殊的GEM算法，它将算法的M步分解成d次条件极大化，每次只改变参数向量的一个分量，其余量不改变。

EM算法另外视角

输入: 观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$

输出: 模型参数$\theta$

1. 选择参数的初值${\theta }^{(0)}$，开始迭代

2. E步：记${\theta }^{(i)}$为第 $i$ 次迭代参数$\theta$的估计值，在第$i+1$次迭代的$E$步，计算
   $$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})=& E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\\ =& \sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array}$$

3. M步
   求使$Q(\theta ,{\theta }^{(i)})$最大化的$\theta$，确定第$i+1$次迭代的参数估计值

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

Q 函数

注意Q函数的定义，可以帮助理解上面E步中的求和表达式

完全数据的对数似然函数$\log P(Y,Z|\theta )$关于给定观测数据$Y$的当前参数${\theta }^{(i)}$下对为观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为Q函数

BMM的EM算法

输入: 观测变量数据$y_1,y_2,\dots ,y_N$，伯努利混合模型

输出: 伯努利混合模型参数

1. 选择参数的初始值开始迭代，$2K$ 个参数
2. E步:
   $${\hat{\gamma }}_{jk}=\frac{{\alpha }_{k}Bern(y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_{k}Bern(y_j|{\theta }_k)}=\frac{{\alpha }_k{\mu }_k^{y_j}(1-{\mu }_k{)}^{1-y_j}}{{\sum }_{k=1}^K{\alpha }_k{\mu }_k^{y_j}(1-{\mu }_k{)}^{1-y_j}},j=1,2,\dots ,N;k=1,2,\dots ,K$$
3. M步:
   $$\begin{gathered} {\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}} \\ {\hat{\alpha }}_k=\frac{n_k}{N} \end{gathered}$$

目标函数L

$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Y,Z|\theta )=\log (\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))$$

目标函数是不完全数据的对数似然

EM算法导出

书中这部分内容回答为什么EM算法能近似实现对观测数据的极大似然估计？
$$\begin{array}{rl}L(\theta )-L({\theta }^{(i)})& =\log \left(\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}\right)-\log P(Y|{\theta }^{(i)})\\ & \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)})\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y|{\theta }^{(i)})\\ & =\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\end{array}$$

以上用于推导迭代过程中两次$L$会变大， 这里面红色部分是后加的方便理解前后两步之间的推导。绿色部分是为了构造期望， 进而应用琴声不等式

高斯混合模型

混合模型，有多种，高斯混合模型是最常用的

高斯混合模型(Gaussian Mixture Model)是具有如下概率分布的模型:
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
其中，${\alpha }_k$是系数，${\alpha }_k\ge 0$，$\sum _{k=1}^K{\alpha }_k=1$, $\varphi (y|{\theta }_k)$ 是高斯分布密度，${\theta }_k=(\mu ,{\sigma }^2)$
$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)$$
上式表示第k个分模型

以上, 注意几点：

1. GMM的描述是概率分布，形式上可以看成是加权求和，有啥用？

2. 加权求和的权重$\alpha$满足${\sum }_{k=1}^K{\alpha }_k=1$的约束

3. 求和符号中除去权重的部分，是高斯分布密度(PDF)。高斯混合模型是一种$\sum (权重\times 分布密度)=分布$的表达
   高斯混合模型的参数估计是EM算法的一个重要应用，隐马尔科夫模型的非监督学习也是EM算法的一个重要应用

4. 书中描述的是一维的高斯混合模型，d维的形式如下[^2]，被称作多元正态分布，也叫多元高斯分布
   $$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{(2\pi {)}^d|\Sigma |}}\exp \left(-\frac{(y-{\mu }_k{)}^T{\Sigma }^{-1}(y-{\mu }_k)}{2}\right)$$
   其中，协方差矩阵$\Sigma \in {\mathbb{R}}^{n\times n}$

5. 另外，关于高斯模型的混合，还有另外一种混合的方式，沿着时间轴的方向做混合。可以理解为滤波器，典型的算法就是Kalman Filter，对应了时域与频域之间的关系，两个高斯的混合依然是高斯，混合的方法是卷积，而不是简单的加法，考虑到的是概率密度的混合，也是一种线性的加权

GMM的EM算法

问题描述:

已知观测数据$y_1,y_2,\dots ,y_N$，由高斯混合模型生成
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
其中， $\theta =({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_K;{\theta }_1,{\theta }_2,\dots ,{\theta }_K)$

补充下，不完全数据的似然函数应该是
$$\begin{array}{rl}P(y|\theta )=& \prod _{j=1}^{N}P(y_j|\theta )\\ =& \prod _{j=1}^N\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\end{array}$$
使用EM算法估计GMM的参数$\theta$

1. 明确隐变量, 初值

• 观测数据$y_j,j=1,2,\dots ,N$这样产生, 是已知的:

  1. 依概率${\alpha }_k$选择第$k$个高斯分布分模型$\varphi (y|{\theta }_k)$;

  2. 依第$k$个分模型的概率分布$\varphi (y|{\theta }_k)$生成观测数据$y_j$

  3. 反映观测数据$y_j$来自第$k$个分模型的数据是未知的, $k=1,2,\dots ,K$ 以**隐变量${\gamma }_{jk}$**表示
     注意这里${\gamma }_{jk}$的维度$(j\times k)$
     γ j k = { 1 , 第 j 个观测来自第 k 个分模型 0 , 否则 j = 1 , 2 , … , N ; k = 1 , 2 , … , K ; γ j k ∈ { 0 , 1 } \gamma_{jk}= \begin{cases} 1, &第j个观测来自第k个分模型\\ 0, &否则 \end{cases}\\ j=1,2,\dots,N; k=1,2,\dots,K; \gamma_{jk}\in\{0,1\} γjk​={1,0,​第j个观测来自第k个分模型否则​j=1,2,…,N;k=1,2,…,K;γjk​∈{0,1}
     注意, 以上说明有几个假设:

  4. 隐变量和观测变量的数据对应, 每个观测数据，对应了一个隐变量，${\gamma }_{jk}$是一种one-hot的形式。

  5. 具体的单一观测数据是混合模型中的某一个模型产生的

• 完全数据为$(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jK},k=1,2,\dots ,N)$

• 完全数据似然函数
  $$\begin{array}{rl}P(y,\gamma |\theta )=& \prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jK}|\theta )\\ =& \prod _{k=1}^K\prod _{j=1}^N{\left[{\alpha }_k\varphi (y_j|{\theta }_k)\right]}^{{\gamma }_{jk}}\\ =& \prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\varphi (y_j|{\theta }_k)\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{(y_j-{\mu }_k{)}^2}{2{\sigma }^2}\right)\right]}^{{\gamma }_{jk}}\end{array}$$
  其中$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$

• 完全数据对数似然函数
  $$\log P(y,\gamma |\theta )=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2}(y_j-{\mu }_k{)}^2\right]\right\}$$

2. E步,确定Q函数

把$Q$ 函数表示成参数形式
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})=& E[\log P(y,\gamma |\theta )|y,{\theta }^{(i)}]\\ =& E\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}]\\ =& E\sum _{k=1}^K\sum _{j=1}^N{\gamma }_{jk}\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}]\\ =& \sum _{k=1}^K\sum _{j=1}^N(E{\gamma }_{jk})\log {\alpha }_k+\sum _{j=1}^N(E{\gamma }_{jk})[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}]\end{array}$$
$$\begin{array}{rl}{\hat{\gamma }}_{jk}=& E({\gamma }_{jk}|y,\theta )=P({\gamma }_{jk}=1|y,\theta )\\ =& \frac{P({\gamma }_{jk}=1,y_j|\theta )}{{\sum }_{k=1}^{K}P({\gamma }_{jk}=1,y_j|\theta )}\\ =& \frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{k=1}^{K}P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}\\ =& \frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)}\end{array}$$

这部分内容就是搬运了书上的公式，有几点说明:

1. 注意这里$E({\gamma }_{jk}|y,\theta )$，记为${\hat{\gamma }}_{jk}$， 对应了E步求的期望中的一部分。
2. 对应理解一下上面公式中的红色，蓝色和绿色部分，以及${\hat{\gamma }}_{jk}$中红色和绿色的对应关系
3. 这里用到了$n_k={\sum }_{j=1}^N{\gamma }_{jk}$
4. ${\hat{\gamma }}_{jk}$为分模型$k$对观测数据$y_j$的响应度。这里，紫色标记的第一行参考伯努利分布的期望。

$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]$$
其中$i$表示第$i$步迭代

1. 写出$Q$ 函数在推导的时候有用，但是在程序计算的时候，E步需要计算的就是${\hat{\gamma }}_{jk}$，M步用到了这个结果。其实抄公式没有什么意义，主要是能放慢看公式的速度。和图表一样，公式简洁的表达了很多信息，公式中也许更能体会到数学之美。

3. M步

求函数$Q(\theta ,{\theta }^{(i)})$对$\theta$的极大值，分别求$\sigma ,\mu ,\alpha$
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

• $\arg \max$ 就是求Q的极值对应的参数$\theta$，如说是离散的，遍历所有值，最大查找，如果是连续的，偏导为零求极值。
• $\frac{\partial Q}{\partial {\mu }_k}=0,\frac{\partial Q}{\partial {\sigma }^2}=0$ 得到${\hat{\mu }}_k,{\hat{\sigma }}_k^2$
• ${\sum }_{k=1}^K{\alpha }_k=1,\frac{\partial Q}{\partial {\alpha }_k}=0$ 得到${\alpha }_k$

4. 停止条件

重复以上计算，直到对数似然函数值不再有明显的变化为止。

如何应用

GMM在聚类中的应用

使用EM算法估计了GMM的参数之后，有新的数据点，怎么计算样本的类别的？

在机器学习[^5]中有一些关于聚类的表述，摘录这里：

Gaussian mixture modeling is among the popular clustering algorithms. The main assumption is that the points, which belong to the same cluster, are distributed according to the same Gaussian distribution(this is how similarity is defined in this case), of unknown mean and covariance matrix.

Each mixture component defines a different cluster. Thus, the goal is to run the EM algorithm over the available data points to provide, after convergence, the posterior probabilities $P(k|x_n),k=1,2,...,K,n=1,2,...,N$, where each k corresponds to a cluster. Then each point is assigned to cluster $k$ according to the rule.

Assign $x_n$to cluster $k=\arg maxP(i|x_n),i=1,2,...,K$

可以参考下scikit-learn的具体实现，就是用了argmax，选择概率最大的那个输出。

Kmeans

另外，直觉上看，GMM最直观的想法就是Kmeans，那么:

1. 在Kmeans常见的描述中都有距离的概念，对应在算法9.2 的描述中，该如何理解?
   这里面距离对应了方差，二范数平方
2. 那么又是怎么在每轮刷过距离之后，重新划分样本的分类呢?
   这里对应了响应度，响应度对应了一个$j\times k$的矩阵，记录了每一个$y_j$ 对第$k$个模型的响应度，可以理解为划分了类别

K怎么定

• 手肘法
• Gap Statistics[^1]
• 第二版中，有对应的描述。$P_{267}$，一般的类别变大的时候，平均直径会增加，当类别数超过某个值之后，平均直径不会变化，这个值可以是最优的$k$