[TOC]
- EM算法的引入
- EM算法
- EM算法的导出
- EM算法在非监督学习中的应用
- EM算法的收敛性
- EM算法在高斯混合模型学习中的应用
- 高斯混合模型
- 高斯混合模型参数估计的EM算法
- EM算法的推广
- F函数的极大极大算法
-
这章如果看三硬币有疑问, 可以往后继续看, 看到高斯混合模型. 然后再回头理解三硬币. 有不理解的地方, 可以重新看对应问题的定义, 重新理解各个符号的意义. 这样也许对学习有帮助.
-
EM算法可以用于生成模型的非监督学习, EM算法是个一般方法, 不具有具体模型.
EM算法是一种迭代算法, 用于含有隐变量的概率模型的极大似然估计,或极大后验概率估计.
-
这里面注意体会不同变量的大小以及对应的取值范围.
-
一个$m\times n\times k$的矩阵可能可以划分成$n$个$m\times k$的形式, 这点理解下.
-
这部分推导有很多求和, 注意体会是按照样本做的, 还是按照模型做的
-
如果对PDF, 高斯分布, 边缘概率分布, 协方差矩阵不清楚, 可以在这个章节从GMM的角度扩展阅读下, 一定会有收获.
-
似然和概率的关系可以推广了解, 这章关于概率和似然的符号表示, 可能会有点看不懂, 比如$P_{157}$中的部分表述. 可以参考引用内容1, 概率和似然是同样的形式描述的都是可能性,
是一个两变量的函数, 似然是给定结果, 求参数可能性; 概率是给定参数求结果可能性. Suppose you have a probability model with parameters
. has two names. It can be called the probability of (given ), or the likelihood of (given that was observed). -
学习过程中注意观测数据在每次EM算法中的意义.
-
GMM中注意区分$\alpha_k$和$\gamma_{jk}$的差异, 直觉上都有一种归属的感觉,
是二值函数, 是一种概率的表示. 是one-hot encoding(also: 1-of-K representation) -
GMM这里面实际上还涉及到一个概念叫做凸组合(Convex Combination)2 . 是凸几何领域的一个概念, 点的线性组合, 所有系数都非负且和为1. 点集的凸包等价于该点集的凸组合.
-
无论是三硬币还是GMM, 采样的过程都是如下:
- Sample
- Sample
注意, 这里用到了$\pi$, 在强化学习中, 随机性策略$\pi(x,a)$表示为状态$x$下选择动作$a$的概率.
- Sample
-
关于EM算法的解释 注意这里EM不是模型, 是个一般方法, 不具有具体的模型.
- PRML
所以, EM应用举例子为kmeans也OK. 而且, 西瓜书$P_{165}$上有说, k均值聚类算法就是一个典型的EM算法
- 统计学习方法
-
函数的极大-极大算法
- PRML
-
这个repo里面实现了BMM算法和GMM算法两种混合模型
-
HMM也是Discrete Dynamic Model, 从图模型角度考虑, 可以发现HMM和卡尔曼滤波以及粒子滤波深层之间的联系. 这部分内容在PRML中有讨论.
一般地, 用$Y$表示观测随机变量的数据,
表示隐随机变量的数据. 和$Z$一起称为完全数据(complete-data), 观测数据$Y$又称为不完全数据(incomplete-data)
上面这个概念很重要, Dempster在1977年提出EM算法的时候文章题目就是《Maximum likelihood from incomplete data via the EM algorithm》, 具体看书中本章参考文献1
假设给定观测数据$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)$
上面这部分简单对应一下, 这里说明一下, 你看到下面概率分布和似然函数形式看起来一样. 在概率中,
观测数据$Y$ | 不完全数据$Y$ | ||
不完全数据$Y$ | 概率分布$P(Y | \theta)$ | 似然函数$P(Y |
完全数据 |
|
\theta )$ | 似然函数$P(Y,Z |
观测数据$Y$
有一点要注意下, 这里没有出现$X$, 在9.1.3节中有提到一种理解
- 有时训练数据只有输入没有对应的输出${(x_1,\cdot),(x_2,\cdot),\dots,(x_N,\cdot)}$, 从这样的数据学习模型称为非监督学习问题.
- EM算法可以用于生成模型的非监督学习.
- 生成模型由联合概率分布$P(X,Y)$表示, 可以认为非监督学习训练数据是联合概率分布产生的数据.
为观测数据, 为未观测数据.
有时候, 只观测显变量看不到关系, 就需要把隐变量引进来.
书中用三硬币模型做为引子, 在学习这部分内容的时候, 注意体会观测数据的作用.
书中用例子来介绍EM算法的问题, 并给出了EM算法迭代求解的过程, 具体例子描述见例9.1, 这块如果不懂, 可以跳过, 看完后面高斯混合模型再回来看.
问题的描述过程中有这样一句: 独立的重复$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)$种可能.
所以上面两个观测序列, 可能出自同一个模型. 在这个问题的求解上是没有区别的, 测试案例$test_t91$做了这个说明. 可以参考.
我们通过一段代码来生成这个数据
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{equation} \begin{aligned} 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{aligned} \end{equation} $$ 以上
- 随机变量$y$是观测变量, 表示一次试验观测的结果是1或0
- 随机变量$z$是隐变量, 表示未观测到的掷硬币$A$的结果
-
是模型参数 - 这个模型是以上数据(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\limits_{Z}P(Z|\theta)P(Y|Z,\theta) $$ 注意这里的求和是下面的"+"描述的部分
即 $$ P(Y|\theta)=\prod\limits^{n}_{j=1}[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}] $$ 注意这里连乘是$Y\rightarrow y_j$出来的, 不理解看似然定义.
考虑求模型参数$\theta=(\pi,p,q)$的极大似然估计, 即 $$ \hat \theta = \arg\max\limits_{\theta}\log P(Y|\theta) $$ 这个题目的标准答案实际上也是未知的. 因为可能生成这样的观测的假设空间太大.
EM算法首选参数初值, 记作$\theta^{(0)}=(\pi^{(0)},p^{(0)}, q^{(0)})$, 然后迭代计算参数的估计值.
如果第$i$次迭代的模型参数估计值为$\theta^{(i)}=(\pi^{(i)}, p^{(i)}, q^{(i)})$
那么第$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{cases} \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{cases} $$ 所以, 这步(求$\mu_j$)干了什么, 样本起到了什么作用?
这一步, 通过假设的参数, 计算了不同的样本对假设模型的响应(
以上, 有点绕.
这一步是什么的期望? 书中有写, **观测数据来自硬币$B$的概率, 在二项分布的情况下, 响应度和概率是一个概念. **这个说明, 有助于后面M步公式的理解.
$$
\begin{align}
\pi^{(i+1)} &= \frac{1}{n}\sum_{j=1}^{n}\mu_j^{(i+1)}\
\color{red}
p^{(i+1)} &= \frac{\sum_{j=1}^{n}\mu_j^{(i+1)}y_j}{\sum_{j=1}^{n}\mu_j^{(i+1)}}\
\color{red}
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{align}
$$
上面, 红色部分的公式从观测数据是来自硬币B的概率
这句来理解.
这个例子里面0.5是个合理又牛逼的初值. 迭代收敛的最后结果是(0.5, 0.6, 0.6)
这个结果说明, 如果A是均匀的, 那么一个合理的解就是B, C是同质的. 他们的分布情况和观测的分布一致.
在测试案例$test_e91$中有计算这部分的结果, 注意看, 这种简单的模型其实收敛的很快.
这里面p对应了A =1, B=1, q对应了A=0, C=1
这三个公式可以改写成如下形式:
$$
\begin{align}
\pi^{(i+1)} &= \frac{1}{n}\sum_{j=1}^{n}\mu_j^{(i+1)}\
\color{red}
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)}\
\color{red}
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{align}
$$
到后面讲高斯混合模型的时候, 可以重新审视这里
$$
\begin{aligned}
\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{aligned}
$$
以上对应了包含两个分量的伯努利混合模型, BMM, 包含四个参数, 因为$\alpha_k$满足等式约束, 所以通常会有三个参数, 另外参见习题$9.3$中有提到两个分量的高斯混合模型的五个参数
实际上也是因为等式约束.
bmm.py对伯努利混合模型做了实现, 有几点说明一下:
-
这个表达式对应了伯努利分布的概率密度, 可以表示成矩阵乘法, 尽量不要用for, 效率会差 -
书中$e_{91}$的表达中, 采用了$\pi, p, q$来表示, 注意在题目的说明部分有说明三个符号的含义
-
实际上不怎么抛硬币, 但是01的伯努利分布很多, 在书中算法9.4部分, 有这样一个说明:
当参数$\theta$的维数为$d(d\ge2 )$的时候, 可以采用一种特殊的GEM算法, 它将算法的M步分解成d次条件极大化, 每次只改变参数向量的一个分量,其余量不改变.
输入: 观测变量数据$Y$, 隐变量数据$Z$, 联合分布$P(Y,Z|\theta)$, 条件分布$P(Z|Y,\theta)$
输出: 模型参数$\theta$
选择参数的初值$\theta^{(0)}$, 开始迭代
E步:记$\theta^{(i)}$为第
次迭代参数$\theta$的估计值, 在第$i+1$次迭代的$E$步,计算 $$ \begin{align} Q(\theta, \theta^{(i)}) =& E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]\ =&\sum_Z\color{red}\log P(Y,Z|\theta)\color{green}P(Z|Y, \theta^{(i)}) \end{align} $$ M步 求使$Q(\theta, \theta^{(i)})$最大化的$\theta$,确定第$i+1$次迭代的参数估计值
注意Q函数的定义
完全数据的对数似然函数$\log P(Y, Z|\theta)$关于给定观测数据$Y$的当前参数$\theta^{(i)}$下对为观测数据$Z$的条件概率分布$P(Z|Y,\theta^{(i)})$的期望称为Q函数.
输入: 观测变量数据$y_1, y_2, \dots, y_N$, 伯努利混合模型
输出: 伯努利混合模型参数
- 选择参数的初始值开始迭代,
个参数 - E步:
- M步:
混合模型, 有多种, 高斯混合模型是最常用的.
高斯混合模型(Gaussian Mixture Model)是具有如下概率分布的模型:
$$
P(y|\theta)=\sum\limits^{K}{k=1}\alpha_k\phi(y|\theta_k)
$$
其中, $\alpha_k$是系数, $\alpha_k\ge0$, $\sum\limits^{K}{k=1}\alpha_k=1$,
以上, 注意几点:
-
GMM的描述是概率分布, 形式上可以看成是加权求和
-
加权求和的权重$\alpha$满足$\sum_{k=1}^K\alpha_k=1$的约束
-
求和符号中除去权重的部分, 是高斯分布密度(PDF). 高斯混合模型是一种$\sum(权重\times 分布密度)=分布$的表达 高斯混合模型的参数估计是EM算法的一个重要应用, 隐马尔科夫模型的非监督学习也是EM算法的一个重要应用.
-
书中描述的是一维的高斯混合模型, d维的形式如下3, 被称作多元正态分布, 也叫多元高斯分布 $$ \phi(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 \R^{n\times n}$
这个弄的不咋好看, plate notation
问题描述:
已知观测数据$y_1, y_2, \dots , y_N$, 由高斯混合模型生成
$$
P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k)
$$
其中,
补充下, 不完全数据的似然函数应该是 $$ \begin{align} P(y|\theta)=&\prod_{j=1}^NP(y_j|\theta)\ =&\prod_{j=1}^N\sum_{k=1}^K\alpha_k\phi(y|\theta_k) \end{align} $$ 使用EM算法估计GMM的参数$\theta$
-
观测数据$y_j, j=1,2,\dots,N$这样产生, 是已知的:
-
依概率$\alpha_k$选择第$k$个高斯分布分模型$\phi(y|\theta_k)$;
-
依第$k$个分模型的概率分布$\phi(y|\theta_k)$生成观测数据$y_j$
-
反映观测数据$y_j$来自第$k$个分模型的数据是未知的,
以**隐变量$\gamma_{jk}$**表示 注意这里$\gamma_{jk}$的维度$(j\times k)$ $$ \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} $$ 注意, 以上说明有几个假设: -
隐变量和观测变量的数据对应, 每个观测数据, 对应了一个隐变量,
是一种one-hot的形式. -
具体的单一观测数据是混合模型中的某一个模型产生的
-
-
完全数据为$(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK},k=1,2,\dots,N)$
-
完全数据似然函数 $$ \begin{aligned} P(y,\gamma|\theta)=&\prod_{j=1}^NP(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK}|\theta)\ =&\prod_{k=1}^K\prod_{j=1}^N\left[\alpha_k\phi(y_j|\theta_k)\right]^{\gamma_{jk}}\ =&\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N\left[\phi(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{aligned} $$ 其中$n_k=\sum_{j=1}^N\gamma_{jk}, \sum_{k=1}^Kn_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} $$
把$Q$ 函数表示成参数形式
$$\begin{aligned}Q(\theta,\theta^{(i)})=&E[\log P(y,\gamma|\theta)|y,\theta^{(i)}]\=&\color{green}E\color{black}\left{\sum_{k=1}^K\left{\color{red}n_k\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N\gamma {jk}\color{black}\left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma k-\frac{1}{2\sigma^2(y_j-\mu_k)^2}\right]\right}\right}\=&\color{green}E\color{black}\left{\sum{k=1}^K\left{\color{red}\sum{j=1}^N\gamma_{jk}\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N\gamma {jk}\color{black}\left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma k-\frac{1}{2\sigma^2(y_j-\mu_k)^2}\right]\right}\right}\=&\sum{k=1}^K\left{\color{red}\sum{j=1}^{N}(\color{green}E\color{red}\gamma_{jk})\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N(\color{green}E\color{blue}\gamma _{jk})\color{black}\left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma _k-\frac{1}{2\sigma^2(y_j-\mu_k)^2}\right]\right}\\end{aligned}$$
$$\begin{aligned}\hat \gamma {jk}= &\color{purple}E(\gamma{jk}|y,\theta)=P(\gamma_{jk}=1|y,\theta)\=&\frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^KP(\gamma_{jk}=1,y_j|\theta)}\=&\frac{P(y_j|\color{red}\gamma_{jk}=1,\theta\color{black})\color{green}P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}\=&\frac{\color{green}\alpha_k\color{black}\phi(y_j|\color{red}\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)}\end{aligned}$$
这部分内容就是搬运了书上的公式, 有几点说明:
- 注意这里$E(\gamma_{jk}|y,\theta)$,记为$\hat\gamma_{jk}$, 对应了E步求的期望中的一部分.
- 对应理解一下上面公式中的红色,蓝色和绿色部分, 以及$\hat\gamma_{jk}$中红色和绿色的对应关系
- 这里用到了$n_k=\sum_{j=1}^N\gamma_{jk}$
-
为分模型$k$对观测数据$y_j$的响应度. 这里, 紫色标记的第一行参考伯努利分布的期望.
- 写出$Q$ 函数在推导的时候有用, 但是在程序计算的时候, E步需要计算的就是$\hat\gamma_{jk}$, M步用到了这个结果.其实抄公式没有什么意义,主要是能放慢看公式的速度. 和图表一样, 公式简洁的表达了很多含义, 公式中也许更能体会到数学之美.
求函数$Q(\theta,\theta^{(i)})$对$\theta$的极大值, 分别求$\sigma, \mu, \alpha$ $$ \theta^{(i+1)}=\arg\max_\theta Q(\theta,\theta^{(i)}) $$
-
就是求Q的极值对应的参数$\theta$, 如说是离散的, 遍历所有值, 最大查找, 如果是连续的, 偏导为零求极值. -
得到$\hat\mu_k, \hat \sigma_k^2$ -
得到$\alpha_k$
重复以上计算, 直到对数似然函数值不再有明显的变化为止.
这部分摘要总结了前面的公式.
因为公式比较集中, 方便对比, 注意体会以下两个方面:
- 这几个公式中待求的变量的维度和角标的关系.
- 这里面有求和, 前面提到过, 注意体会每一步刷的是模型, 还是样本
另外, 直觉上看, GMM最直观的想法就是Kmeans, 那么:
- 在Kmeans常见的描述中都有距离的概念, 对应在算法9.2 的描述中, 该如何理解? 这里面距离对应了方差
- 那么又是怎么在每轮刷过距离之后, 重新划分样本的分类呢? 这里对应了响应度, 响应度对应了一个$j \times k$的矩阵, 记录了每一个$y_j$ 对第$k$个模型的响应度, 可以理解为划分了类别.
- 手肘法
- Gap Statistics4
广义期望极大(generalized expectation maximization,
- 关于习题9.3 GMM模型的参数($\alpha _k, \mu k, \sigma^2_k $)应该是$3k$个, 题目9.3中提出两个分量的高斯混合模型的5个参数, 是因为参数$\alpha_k$满足$\sum{k=1}^K\alpha _k=1$