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

文章目录

• • 导读
  • 符号说明
  • 混合模型
  • • 伯努利混合模型(三硬币模型)
    • • 问题描述
      • 三硬币模型的EM算法
      • • 1.初值
        • 2.E步
        • 3.M步
        • 初值影响
        • p,q 含义
    • EM算法另外视角
    • • Q 函数
      • BMM的EM算法
      • 目标函数L
      • EM算法导出
    • 高斯混合模型
    • • GMM的EM算法
      • • 1. 明确隐变量, 初值
        • 2. E步,确定Q函数
        • 3. M步
        • 4. 停止条件
      • 如何应用
      • • GMM在聚类中的应用
    • Kmeans
    • • K怎么定

导读

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

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

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

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

  注意，这里用到了$\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作了两个基本假设，实际上是在说在图模型中，存在哪些边

符号说明

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

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

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

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

混合模型

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

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

问题描述

问题的描述过程中有这样一句：独立的重复$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 np

p = 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 = 0
        print("[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 = 0
        print("[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$, 则观测数据的似然函数为

即
$$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的伯努利分布很多

EM算法另外视角

Q 函数

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

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

BMM的EM算法

目标函数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]中有一些关于聚类的表述，摘录这里：

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

Kmeans

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

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

K怎么定

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