Expectation Maximum Algorithm(EM

作者: To_QT | 来源:发表于2019-08-16 14:41 被阅读2次

1. 预备知识

1.1 凸函数的性质

假设定义在实数域上的函数，对于任意的实数，都有
$f''(x)<0\tag{1.1}$
则函数 $f(x)$ 称为凸函数，反之，为凹函数。（国内外教材对于凸凹函数的定义不同，这里采用国内教材为准）

1.2 Jensen不等式

1.2.1 数学意义

如果函数是凸函数，则有
$f(\lambda x_1+(1-\lambda)x_2) \geq \lambda f(x_1)+(1-\lambda)f(x_2)\tag{1.2}$

如果函数是凹函数，则有
$f(\lambda x_1+(1-\lambda)x_2) \leq \lambda f(x_1)+(1-\lambda)f(x_2)\tag{1.3}$

1.2.2 几何意义

以凹函数为例

当 $\lambda$ 为 $0.5$ 时，左右两边严格相等。

1.2.3 推广

若存在凸函数 $f(x)$ ，其中有
$\begin{align} \sum_{i=1}^{k} \lambda_i=1, \lambda_{i} \geq 0\tag{1.4} \end{align}$
则有
$\begin{align} f(\lambda_1x_1+\lambda_2x_2+...+\lambda_kx_k)\geq\lambda_{1}f(x_1)+\lambda_{2}f(x_2)+...+\lambda_{k}f(x_k)\tag{1.5} \end{align}$
在概率论中，我们知道对于离散变量的期望 $E$ 有
$\begin{align} E(X)=\sum_{i}^{k}p_iX_i \tag{1.6} \end{align}$
类比公式 $(1.4),(1.5)$ 则有
$\begin{align} f(E(X)) \geq E(f(X)) \tag{1.7} \end{align}$

1.3 高斯分布与高斯混合分布

1.3.1 高斯分布

说“高斯分布”可能有点不适应，但是，其实高斯分布（Gaussian distribution）就是正态分布（Normal distribution），也称“常态分布”。
若随机变量 $X$ （总共有 $N$ 个样本， $x_1，x_2…,x_N$ ）服从一个位置参数为 $\mu$ 、尺度参数为 $\sigma$ 的概率分布，且其概率密度函数为
$\begin{align} f(x)=\frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(x-\mu)^2}{2\sigma^2})\tag{1.8} \end{align}$
则称为正态分布。

1.3.2 高斯混合分布

通常在自然界中，随机变量 $X$ 可能有 $K$ 个高斯分布混合组成，那么，记取自不同高斯分布的概率分别为 $\pi_1, \pi_2,... ,\pi_K$ ，第 $k$ 个高斯分布的均值为 $\mu_k$ ，方差为 $\sigma_k$ 。

注意一点：

当随机变量 $X$ 是 $1$ 维数据时，均值 $\mu_k$ 和方差 $\sigma_k$ 均为一个标量（其中 $(k<K)$ ）

当随机变量 $X$ 是 $m$ 维数据时，均值 $\mu_k$ 是一个 $m$ 维的向量，和方差 $\sigma_k$ 是一个 $m*m$ 的协方差矩阵。为了加以区分，将 $\sigma_k$ 记为 $\Sigma_k$ 。

此时，需要估计的参数包括 $K$ 个 $\pi, \mu, \Sigma$ 。
因此，高斯混合模型是指具有如下形式的概率分布模型：
$\begin{align} p(y|\theta)=\sum_{k=1}^{K}a_k\phi(y; \theta_k)\tag{1.9} \end{align}$
其中， $a_k$ 是系数， $a_k \geq 0, \sum_{k=1}^{K}a_k=1$ ，且 $\phi(y|\theta_k)$ 是高斯分布密度 $\phi(y; \mu_k, \Sigma_k)$ 。

1.4 隐变量

不能被直接观察到，但是对系统的状态和能观察到的输出存在影响的一种东西。

比如投硬币，规则如下：甲投掷硬币，乙猜测硬币的正反面。甲投硬币的规则如下：A、B、C三枚硬币，A正则选择B，A反选C，B正为1，B反为0；C正为1，C反为0。那么，假设随机变量 $y$ 表示一次试验的观测结果： $0$ 或 $1$ ，随机变量 $z$ 表示投掷A硬币的结果。如果在实验过程中，乙知道A硬币的结果，则 $z$ 不为隐变量；若乙不知道硬币A的投币结果，则 $z$ 为隐变量。

1.5 极大似然函数的概念

1.5.1 极大似然函数的概念

1.5.2 极大似然函数——以高斯分布为例

若给定一组样本 $X_1，X_2…,X_N$ ，已知它们来自于高斯分布 $N~(\mu, \sigma^2)$ ，试估参数 $\mu, \sigma$ 。

1.5.3 极大似然函数——以高斯混合分布为例

对于高斯混合分布，建立似然函数，分两步走：

因为第 $i$ 个样本属于第 $k$ 个分布的概率为： $\pi_kN(X_i; \mu_k, \Sigma_k)$ ，所以取到第 $i$ 个样本的概率为：
$\begin{align} \sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k)\tag{1.10} \end{align}$
所以， $N$ 个样本的似然函数为
$\begin{align} L_{\pi, \mu, \Sigma}(X)=&\prod_{i=1}^{N} \sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k) \tag{1.11} \end{align}$
对公式 $(1.11)$ 取对数有
$\begin{align} log(L_{\pi, \mu, \Sigma}(X_i))=&log(\prod_{i=1}^{N} \sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k)))\\ =&\sum_{i=1}^{N} log(\sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k)))\\ \tag{1.12} \end{align}$
此时，我们的目标是要找到 $K$ 个高斯分布的 $\mu, \Sigma$ ，使得公式 $(1.12)$ 最大
$\begin{align} log(L_{\pi, \mu, \Sigma}(X_i))=&\sum_{i=1}^{N} log(\sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k)))\\ \tag{1.13} \end{align}$

2. EM算法的直观理解

先假设参数，根据参数计算其所属概率，反过来计算已知概率情况下该参数的后验概率，再通过最大似然估计参数，循环往复，直到收敛。

以知乎男女身高分布的例子为例。随机挑选1000个人，测量他们的身高。在这1000个人中，有男性有女性，身高分别服从 $N_{boy}(\mu_{boy}, \sigma_{boy}^2)$ 和 $N_{girl}(\mu_{girl}, \sigma_{girl}^2)$ 的分布，试着估计 $\mu_{boy}, \sigma_{boy}, \mu_{girl}, \sigma_{girl}$ 。

估计每个样本数据由每个分布组成的比例。对于每个样本 $X$ ，它由第 $k$ 个高斯组份所占的比例为

$\begin{align} \gamma(i,k)=\frac{\pi_kN(X_i;\mu_k,\sigma_k)}{\sum_{k=1}^{K}\pi_kN(X_i;\mu_k,\sigma_k)} \tag{2.1} \end{align}$
首先假设男孩服从的分布为 $N_{boy}(175, 10^2)$ ，女孩服从分布 $N_{girl}(165, 8^2)$ 。将第一个样本 $X_1=198$ 分别带入公式 $(2.2-1)(2.2-2)$ 中，有：
$\begin{align} f_{boy}(198)=\frac{1}{\sqrt{2\pi}*10}e^{-\frac{(198-175)^2}{2*10^2}}\tag{2.2-1} \end{align}$
$\begin{align} f_{girl}(198)=\frac{1}{\sqrt{2\pi}*8}e^{-\frac{(198-165)^2}{2*8^2}}\tag{2.2-2} \end{align}$

其次，假设选择男生分布的概率为 $\pi_{boy}=\frac{1}{2}$ ，女性分布的概率为 $\pi_{girl}=\frac{1}{2}$ ，带入公式 $(2.1)$ 可以写成

$\begin{align} \gamma(1,boy)=&\frac{\pi_{boy}N(X_1; \mu_{boy},\sigma_{boy})}{\sum_{k=1}^{K}\pi_kN(X_1;\mu_k,\sigma_k)} \\\\ =&\frac{\frac{1}{2}*f_{boy}(198)}{\frac{1}{2}*f_{boy}(198)+\frac{1}{2}*f_{girl}(198)}\tag{2.3-1} \end{align}$
$\begin{align} \gamma(1,girl)=&\frac{\pi_{girl}N(X_1; \mu_{girl},\sigma_{girl})}{\sum_{k=1}^{K}\pi_kN(X_1;\mu_k,\sigma_k)} \\\\ =&\frac{\frac{1}{2}*f_{girl}(198)}{\frac{1}{2}*f_{boy}(198)+\frac{1}{2}*f_{girl}(198)}\tag{2.3-2} \end{align}$
此时，男孩样本 $1$ 上的变化情况为 $198*\gamma(1,boy)$ ，女孩样本 $1$ 上的变化情况为 $198*\gamma(1,girl)$ 。

对所有样本均重复上述操作。

可得到 $\gamma(i,boy)$ 、 $\gamma(i,girl)$ ，其中 $i=1,2,...,N$ 。第一轮参数估计的结果为：
$\begin{align} \mu_{boy}=&\frac{1}{\sum_{i=1}^{N}\gamma(i,boy)}\sum_{i=1}^{N}X_i*\gamma(i,boy) \tag{2.4-1} \\ \mu_{girl}=&\frac{1}{\sum_{i=1}^{N}\gamma(i,girl)}\sum_{i=1}^{N}X_i*\gamma(i,girl) \tag{2.4-2} \\ \sigma_{boy}=&\frac{1}{\sum_{i=1}^{N}\gamma(i,boy)}\sum_{i=1}^{N}\gamma(i,boy)*(X_i-\mu_{boy})^2 \tag{2.4-3} \\ \sigma_{girl}=&\frac{1}{\sum_{i=1}^{N}\gamma(i,girl)}\sum_{i=1}^{N}\gamma(i,girl)* (X_i-\mu_{girl})^2 \tag{2.4-4} \\ \pi_{boy}=&\frac{\sum_{i=1}^{N}*\gamma(i,boy)}{N} \tag{2.4-5} \\ \pi_{girl}=&\frac{\sum_{i=1}^{N}*\gamma(i,girl)}{N} \tag{2.4-6} \\ \end{align}$

根据公式 $(2.4-1,2.4-2,2.4-3,2.4-4, 2.4-5,2.4-6)$ 所得，重复上述步骤，直到迭代次数完毕为止。

3. EM算法的科学推理

3.1 EM算法问题提出

对于不含隐变量的数据而言，多元高斯分布的似然函数可以写成
$\begin{align} log(L_{\pi, \mu, \Sigma}(X_i))=&\sum_{i=1}^{N} log(\sum_{k=1}^{K} \pi_kN(X_i; \mu_k, \Sigma_k)))\\ \tag{1.13} \end{align}$
一般情况下，我们会对 $\mu_k,\Sigma_k$ 求偏导，但是， $(3.3)$ 中的函数跟隐变量没有半毛钱关系，但是实际生活中，我们却不得不考虑隐变量，因此，需要稍微变动一下。

3.2 加入了隐变量之后的似然函数

其实前面都是瞎比比，没有卵用，这里才是主题。。。

首先假设，训练集由 $K$ 个高斯混合分布组成。在训练集 $\boldsymbol{X}=[\boldsymbol{X_1}, \boldsymbol{X_2}, ..., \boldsymbol{X_N}]^T$ （对应上例中的身高）中，给定每个样本具有 $d$ 个特征值。除此之外，还具有未观测数据 $\boldsymbol{z}=[z_1, z_2, ..., z_N]^T$ 。所以，未观测数据所构成的分布为 $Q_i(z_i=K)$ ，可能的取值（随机变量）为： $z_i \in \left \{1, 2, ..., K\right \}$ ，例如，在第 $i$ 个样本中， $\boldsymbol{X_i}|(Z_i=1)\sim N_d(\boldsymbol{\mu_1}, \Sigma_1)$ ； $\boldsymbol{X_i}|(Z_i=2)\sim N_d(\boldsymbol{\mu_2}, \Sigma_2)$ ；...； $\boldsymbol{X_i}|(Z_i=K)\sim N_d(\boldsymbol{\mu_K}, \Sigma_K)$ 。表示生成样本 $\boldsymbol{X_i}$ 的高斯混合成分，取值未知。（对应上例中的性别）。

需要解决的问题：

隐变量与高斯混合分布的关系是怎样的？

隐变量与高斯混合分布是一一对应的，隐变量反映了观测数据来自K个高斯混合分布中的哪一个分布。在身高例子中， $z_i=boy$ 说明，第 $i$ 个变量来自于“男”这个高斯分布。

隐变量怎么与公式 $(1.13)$ 结合?

隐变量与观测变量构成联合概率，因此，第 $i$ 个样本的隐变量与观测变量构成的分布为 $p(\boldsymbol{X_i}, z_i; \theta, \Sigma_i)$

3.2.1 明确隐变量，写出似然函数

$\begin{align} l(\theta)=&\sum_{i=1}^{N} log(\sum_{\boldsymbol{Q_i(Z)}}\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta))\\ \tag{2.5} \end{align}$

一般情况下，由于 $z_i$ 是隐随机变量，不便于直接寻找其参数估计，因此，采用的策略是：计算 $l(\theta)$ 的下界，求其下界的最大值，不断重复这个过程，直到找到极值为止。也就是说，其终止条件为当下界函数的极大值与 $l(\theta)$ 在该点的取值相等时，则说明该点为 $l(\theta)$ 的极大值。

首先 $Q$ 是隐随机变量 $z_i$ 的某一个分布，且 $Q(z_i)>0$ ，则
$\begin{align} l(\theta)=&\sum_{i=1}^{N} log(\sum_{\boldsymbol{Q_i(Z)}}\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta))\\ =&\sum_{i=1}^{N} log(\sum_{\boldsymbol{Q_i(Z)}}Q_i (z_i)\frac{\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta)}{Q(z_i)})\\ \tag{2.6} \end{align}$
由于 $log$ 函数是凸函数，则根据公式 $(1.5)$ ，公式 $(2.6)$ 可以改为：
$\begin{align} l(\theta)=&\sum_{i=1}^{N} log(\sum_{\boldsymbol{Q_i(Z)}}\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta))\\ =&\sum_{i=1}^{N} log(\sum_{\boldsymbol{Q_i(Z)}}Q_i(z_i)\frac{\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta)}{Q_i(z_i)})\\ \geq& \sum_{i=1}^{N} \sum_{\boldsymbol{Q_i(Z)}}Q_i(z_i) log(\frac{\pi_{z_i}p(\boldsymbol{X_i},z_i;\theta)}{Q_i(z_i)}) \tag{2.7} \end{align}$

2.2.2 如何构建并求解下界函数？

E-step：

2.1.1节中提到，当下界函数的极大值与 $l(\theta)$ 在该点的取值相等时，则说明该点为 $l(\theta)$ 的极大值。那么，什么时候会出现这种情况呢？
如果公式 $(1.5)$ 中的变量 $x$ 是常数呢，那么，此时等号也成立。此时有
$\frac{p(\boldsymbol{X_i},z_i;\theta)}{Q_i(z_i)}=C,\ C是一个常数。\tag{2.8}$
则公式 $(2.8)$ 可以改写成
$p(\boldsymbol{X_i},z_i;\theta)=Q_i(z_i)*C\tag{2.9}$
两边同时累加可得
$\begin{align} \sum_{\boldsymbol{Q_i(Z)}}p(\boldsymbol{X_i},z_i;\theta)=&\sum_{\boldsymbol{Q_i(Z)}}Q_i(z_i)*C\\ \sum_{\boldsymbol{Q_i(Z)}}p(\boldsymbol{X_i},z_i;\theta)=&C \tag{2.10} \end{align}$
结合公式 $(2.8)$ 与公式 $2.10$ 可知：
$\begin{align} Q_i(z_i)=&\frac{p(\boldsymbol{X_i},z_i;\theta)}{\sum_{\boldsymbol{Q_i(Z)}}p(\boldsymbol{X_i},z_i;\theta)} \\ =&\frac{p(\boldsymbol{X_i},z_i;\theta)}{p(\boldsymbol{X_i};\theta)} \tag{2.11} \end{align}$
所以，若 $Q_i(z_i)$ 为给定第 $i$ 个样本 $\boldsymbol{X_i}$ 的条件下，随机变量 $\boldsymbol{Z}=z_i$ 的条件概率时，等号成立。(对于身高198的人来说，若求的是男生的参数估计时，需要计算在身高198的条件下，性别是男的概率。

M-step：将参数带入：

$\begin{align} \sum_{i=1}^{N}\sum_{\boldsymbol{Q_i(Z)}}Q_i(z_i)log(\frac{p(\boldsymbol{X_i},z_i;\theta)}{Q_i(z_i)}) \\ \sum_{i=1}^{N} \sum_{\boldsymbol{Q_i(Z)}}Q_i(z_i) log(\frac{p(\boldsymbol{X_i}|z_i;\theta)p(z_i;\theta)}{Q_i(z_i)}) \\ \sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\frac{p(\boldsymbol{X_i}|z_i=j;\theta)p(z_i=j;\theta)}{Q_i(z_i=j)})\tag{2.12} \end{align}$
那么在高斯混合分布中则有：
$\begin{align} \sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\frac{\frac{1}{{(2\pi)^{n/2}\boldsymbol{|\Sigma_j|^{1/2}}}}exp(-\frac{1}{2}(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{X_i}-\boldsymbol{\mu_j})\pi_j}{Q_i(z_i=j)})\tag{2.13} \end{align}$
公式 $(2.13)$ 中:

向量 $\boldsymbol{X_i}$ 代表第 $i$ 个样本的特征数据；
向量 $\boldsymbol{\mu_i}$ 代表第 $j$ 个分布的均值；
矩阵 $\boldsymbol{\Sigma_i}$ 代表第 $j$ 个分布的协方差矩阵；
标量 $z_j$ 代表第 $i$ 个样本中隐随机变量的取值（男or女）
标量 $\pi_j$ 代表第 $j$ 个分布 $z_i$ 所对应的先验概率

2.2.3 计算参数 $\mu_j$ 、 $\Sigma_j$ 、 $\pi_j$

对参数 $\mu_j$ 、 $\Sigma_j$ 分别求导，矩阵求导参考。

$\partial{\boldsymbol{X^TAX}}/\partial{\boldsymbol{X}}=2\boldsymbol{AX}$ ， $\boldsymbol{A}$ 为对称矩阵。

计算参数 $\mu_j$ 。
$\begin{align} &= \bigtriangledown_{\mu_j}\sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\frac{\frac{1}{{(2\pi)^{n/2}\boldsymbol{|\Sigma_j|^{1/2}}}}exp(-\frac{1}{2}(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{X_i}-\boldsymbol{\mu_j})\pi_j}{Q_i(z_i=j)})\\ &= \bigtriangledown_{\mu_j} \sum_{i=1}^{N} \sum_{j=1}^{K}-\frac{1}{2}Q_i(z_i=j)(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{X_i}-\boldsymbol{\mu_j})\\ &=\sum_{i=1}^{N}Q_i(z_i=j)\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{\mu_j}-\boldsymbol{X_i}) \end{align}$
令求导结果为 $0$ ，则有
$\boldsymbol{\mu_j}=\frac{\sum_{i=1}^{N}Q_i(z_i=j)(\boldsymbol{X_i})}{\sum_{i=1}^{N}Q_i(z_i=j)}\tag{2.14-1}$
计算参数 $\Sigma_j$ 。
$\begin{align} &=\bigtriangledown_{\Sigma_j}\sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\frac{\frac{1}{{(2\pi)^{n/2}\boldsymbol{|\Sigma_j|^{1/2}}}}exp(-\frac{1}{2}(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{X_i}-\boldsymbol{\mu_j})\pi_j}{Q_i(z_i=j)})\\ &= \bigtriangledown_{\Sigma_j} \sum_{i=1}^{N} \sum_{j=1}^{K}-\frac{1}{2}Q_i(z_i=j)(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T\boldsymbol{\Sigma_j}^{-1}(\boldsymbol{X_i}-\boldsymbol{\mu_j})\\ {\Sigma_j} &= \frac{\sum_{i=1}^{N}Q_i(z_i=j)(\boldsymbol{X_i}-\boldsymbol{\mu_j})(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T}{\sum_{i=1}^{N}Q_i(z_i=j)} \tag{2.14-2} \end{align}$
计算参数 $\pi_j$ 。
在 $\mu_j$ 和 $\Sigma_j$ 已知的情况下，且 $\pi_j \geq 0$ ，有 $\sum_{j=1}^{K}\pi_j=1$ 。因此，目标函数 $(2.13)$ 可化简为：
$\begin{align} \sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\pi_j) \tag{2.15} \end{align}$
引入拉格朗日算子：
$\begin{align} l(\pi_j)=\sum_{i=1}^{N} \sum_{j=1}^{K}Q_i(z_i=j) log(\pi_j)+\lambda(\sum_{j=1}^{K}\pi_j-1) \tag{2.16} \end{align}$
对 $\pi_j$ 求偏导，则有：
$\begin{align} \frac{\partial{l(\pi_j)}}{\partial{\pi_j}}=\sum_{i=1}^{N}Q_i(z_i=j) \frac{1}{\pi_j}+\lambda \tag{2.17} \end{align}$
令公式 $(2.17)=0$ ，则有：
$\begin{align} 0&=\sum_{i=1}^{N}Q_i(z_i=j) \frac{1}{\pi_j}+\lambda\\ 0&=\sum_{i=1}^{N}Q_i(z_i=j) +\lambda\pi_j\\ \pi_j&=-\frac{1}{\lambda}\sum_{i=1}^{N}Q_i(z_i=j)\tag{2.17-1}\\ 0&=\sum_{j=1}^{K}\sum_{i=1}^{N}Q_i(z_i=j)+\lambda\sum_{j=1}^{K}\pi_j\\ 0&=\sum_{i=1}^{N}1+\lambda\\ \lambda &= -N \tag{2.17-2} \end{align}$
所以
$\begin{align} \pi_j&=\frac{1}{M}\sum_{i=1}^{N}Q_i(z_i=j)\\ \tag{2.18} \end{align}$
综上所述：
$\begin{align} \boldsymbol{\mu_j}&=\frac{\sum_{i=1}^{N}Q_i(z_i=j)(\boldsymbol{X_i})}{\sum_{i=1}^{N}Q_i(z_i=j)}\\ {\Sigma_j} &= \frac{\sum_{i=1}^{N}Q_i(z_i=j)(\boldsymbol{X_i}-\boldsymbol{\mu_j})(\boldsymbol{X_i}-\boldsymbol{\mu_j})^T}{\sum_{i=1}^{N}Q_i(z_i=j)}\\ \pi_j&=\frac{1}{M}\sum_{i=1}^{N}Q_i(z_i=j)\\ \end{align}$