第九章-EM算法

从第九章开始，学习总结的东西有所不同了，第2-8章是分类问题，都属于监督学习，第9章EM算法是非监督学习。本文主要是总结EM算法的应用以及处理问题的过程和原理推导。

EM算法

EM算法(期望极大算法 expectation maximization algorithm)是一种迭代算法。当我们面对概率模型的时候，既有观测变量，又含有隐变量或者潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接使用极大似然估计法或者贝叶斯估计模型估计参数，但是，当模型含有隐变量的时候，就不能简单地这样估计，此时，在1977年，Dempster等人总结提出EM算法：E步：求期望(expectation);M步：求极大值(maximization)。
$$\begin{gathered} 输入：观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|\theta )，条件分布P(Z|Y,\theta )。 \\ 输出：模型参数\theta 。 \\ (1)选择参数的初值{\theta }^{(0)}，开始迭代。 \\ (2)\ast \ast E步：\ast \ast 记{\theta }^{(i)}为第i次迭代参数\theta 的估计值，在第i+1次迭代的E步， \\ 计算\begin{array}{rlr}Q(\theta ,{\theta }^{(i)})=& E_Z[\ln P(Y,Z|\theta )|Y,{\theta }^{(i)}]=& \sum _Z\ln P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\end{array} \\ 这里，P(Z|Y,{\theta }^{(i)})是在给定观测数据Y和当前的参数估计{\theta }^{(i)}下隐变量数据Z的条件概率分布。 \\ (3)\ast \ast M步：\ast \ast 求使得Q(\theta ,{\theta }^{(i)})极大化的\theta ，确定第i+1次迭代的参数估计值{\theta }^{(i+1)}{\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)}) \\ (4)重复第(2)步和第(3)步，直到收敛（收敛条件：{\theta }^{(i)}和{\theta }^{(i+1)}很接近， \\ 或者是Q({\theta }^{(i+1)},{\theta }^{(i)})和Q({\theta }^{(i)},{\theta }^{(i-1)})很接近）。函数Q(\theta ,{\theta }^{(i)})是EM算法的核心，称为Q函数。 \end{gathered}$$

推导过程

上述阐述了EM算法，可是为什么EM算法能近似实现对观测数据的极大似然估计呢？下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法，从而了解EM算法的作用。

在推导过程中用到的公式：

$$\begin{gathered} Jenson不等式：f(\sum _i{\alpha }_i{x}_i)⩾\sum _i{\alpha }_{i}f(x_i)其中函数f是凸函数， \\ 那么对数函数也是凸函数，\sum _i{\alpha }_i=1，{\alpha }_i表示权值，0⩽{\alpha }_i⩽1 \end{gathered}$$

$$\begin{gathered} 首先有一个需要观测的向量\theta ，观测数据Y=(y_1,y_2,\cdots ,y_N)，隐变量Z=(z_1,z_2,\cdots ,z_N)， \\ 当求解\theta 时，似然函数为\begin{array}{r}L(\theta )=\ln P(Y|\theta )=\ln \sum _{Z}P(Y,Z|\theta )=\ln (\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))\end{array} \\ 假设在第i次迭代后\theta 的估计值为{\theta }^{(i)}，希望新估计值\theta 能使L(\theta )增加，即L(\theta )>L({\theta }^{(i)})，则可计算两者的差： \\ L(\theta )-L({\theta }^{(i)})=\ln (\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))-\ln P(Y|{\theta }^{(i)}) \\ 一般来说，对\ln P_1{P}_2\cdots P_N比较好处理，但是如果是\ln \sum P_1{P}_2就不好处理， \\ 为了将\sum 求和符号去掉，用Jenson不等式进行缩放处理。 \\ 对于上述形式，对Z求和，要如何凑出来一个具有Jenson不等式中的{\alpha }_i呢？ \\ 很容易想到，关于Z的密度函数，该密度函数取值求和为1，需要构造一个Z的概率分布。 \\ \begin{array}{r}L(\theta )-L({\theta }^{(i)})=\ln (\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))-\ln P(Y|{\theta }^{(i)})\\ =\ln (\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})})-\ln P(Y|{\theta }^{(i)})\end{array} \\ 由Jesson不等式，\ln (\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})})⩾\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})} \\ \because \ln P(Y|{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln P(Y|{\theta }^{(i)}) \\ \begin{array}{r}\therefore L(\theta )-L({\theta }^{(i)})⩾\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln P(Y|{\theta }^{(i)})\\ =\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\end{array} \\ 令B(\theta ,{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})} \\ \therefore L(\theta )⩾B(\theta ,{\theta }^{(i)})， \\ 也就是说B(\theta ,{\theta }^{(i)})是L(\theta )的一个下界，要最大化L(\theta )，即最大化B(\theta ,{\theta }^{(i)})。 \\ \begin{array}{r}\therefore {\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}B(\theta ,{\theta }^{(i)})=\underset{\theta }{\mathrm{argmax}}(\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln P(Z|\theta )P(Y|Z,\theta ))\\ =\underset{\theta }{\mathrm{argmax}}(\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln P(Y,Z|\theta ))\end{array} \\ \because Q(\theta ,{\theta }^{(i)})=\sum _Z\ln P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)}) \\ \therefore {\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}(Q(\theta ,{\theta }^{(i)})) \\ 等价于EM算法的M步，E步等价于求\sum _{Z}P(Z|Y,{\theta }^{(i)})\ln P(Y,Z|\theta )， \\ 以上就得到了EM算法，通过不断求解下界的极大化逼近求解对数似然函数极大化。 \end{gathered}$$

EM算法在高斯混合模型学习中的应用

高斯混合模型

$$\begin{gathered} 高斯混合模型是指具有如下形式的概率分布模型： \\ P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k) \\ 其中，{\alpha }_k是系数，{\alpha }_k⩾0,\sum _{k=1}^K{\alpha }_k=1，\varphi (y|{\theta }_k)是高斯分布密度， \\ {\theta }_k=({\mu }_k,{\sigma }_k^2)，\varphi (y|\theta )=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)称为第k个分模型。 \\ 首先介绍高斯混合模型，只考虑简单的一维随机变量y，高斯分布就是正态分布， \\ y\sim N(\mu ,{\sigma }^2)，给定y的观测值，就可以很容易求出\mu 和{\sigma }^2，但是目前y不是来自高斯分布， \\ 而是有一定概率的来自于两个不同的高斯分布N({\mu }_1,{\sigma }_1^2)和N({\mu }_2,{\sigma }_2^2)，这个就是两个高斯分布的混合， \\ 并不知道y来自于哪一个高斯分布，这里涉及到了隐变量。对于包含隐变量的参数估计，对此可以做以下处理。 \\ 用一个向量\gamma 表示z，如果z=1，则\gamma =(1,0,0,\cdots ,0)， \\ 如果z=2，则\gamma =(0,1,0,\cdots ,0)，这个相当于One-hot， \\ 也就是说z是第i个高斯分布，在\gamma 的第i个分量为1，其他分量都为0。 \end{gathered}$$

推导过程

明确隐变量，写出完全数据的对数似然函数

根 据 E M 算 法 ， 存 在 一 个 隐 变 量 γ ， γ 表 示 当 前 的 y 来 自 的 高 斯 分 布 ， 对 于 第 一 个 观 测 值 ， 有 γ 1 = ( γ 11 , γ 12 , ⋯   , γ 1 K ) ， 其 中 根 据 书 中 的 γ j k 的 定 义 ： γ j k = { 1 , 第 j 个 观 测 来 自 第 k 个 分 模 型 0 , 否 则 j = 1 , 2 , ⋯   , N ; k = 1 , 2 , ⋯   , K 以 上 是 随 机 变 量 的 分 布 ， 取 第 1 个 值 的 概 率 为 α 1 ， 取 第 2 个 值 的 概 率 为 α 2 ， … … ， 取 第 K 个 值 的 概 率 为 α K ， 一 旦 知 道 γ 1 的 值 ， 就 知 道 从 第 几 个 高 斯 分 布 中 抽 取 y 1 。 p ( γ 1 , y 1 ∣ θ ) = p ( γ 1 ∣ θ ) ⋅ p ( y 1 ∣ γ 1 , θ ) = α γ 11 1 ⋅ α γ 12 2 ⋯ α γ 1 K K ϕ ( y 1 ∣ θ 1 ) γ 11 ϕ ( y 2 ∣ θ 2 ) γ 12 ⋯ ϕ ( y 1 ∣ θ K ) γ 1 K = ∏ k = 1 K [ α k ϕ ( y 1 ∣ θ k ) ] γ 1 k 这 个 是 第 1 个 样 本 点 完 全 数 据 的 密 度 函 数 。 在 极 大 化 似 然 估 计 中 是 极 大 化 似 然 函 数 ， 这 需 要 所 有 样 本 点 的 联 合 分 布 ， 对 于 所 有 的 样 本 点 ， 概 率 密 度 函 数 为 P ( y , γ ∣ θ ) = ∏ j = 1 N ∏ k = 1 K [ α k ϕ ( y j ∣ θ k ) ] γ j k ∵ ∏ j = 1 N ∏ k = 1 K α k γ j k = ∏ k = 1 K α k ∑ j = 1 N γ j k ， ∑ j = 1 N γ j k 表 示 在 N 个 样 本 点 中 ， 一 共 有 多 少 个 是 来 自 第 k 个 高 斯 分 布 的 ， 将 该 数 量 记 为 n k = ∑ j = 1 N γ j k ， n 1 + n 2 + ⋯ + n K = N ∴ ∏ k = 1 K α k ∑ j = 1 N γ j k = ∏ k = 1 K α k n k ∴ ∏ j = 1 N ∏ k = 1 K α k γ j k = ∏ k = 1 K α k n k ∴ P ( y , γ ∣ θ ) = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y i ∣ θ k ) ] γ j k = ∏ k = 1 K α k n k ⋅ ∏ j = 1 N [ 1 2 π σ k exp ⁡ ( − ( y j − μ k ) 2 2 σ k 2 ) ] γ j k ∴ ln ⁡ P ( y , γ ∣ θ ) = ∑ k = 1 K { n k ln ⁡ α k + ∑ j = 1 N γ j k [ ln ⁡ ( 1 2 π ) − ln ⁡ σ k − 1 2 σ k 2 ( y j − μ k ) 2 ] } 根据EM算法，存在一个隐变量\gamma，\gamma表示当前的y来自的高斯分布，对于第一个观测值，\\有\gamma_1=(\gamma_{11},\gamma_{12},\cdots,\gamma_{1K})，其中根据书中的\gamma_{jk}的定义：\gamma_{jk} = \left\{ \begin{aligned} 1, & 第j个观测来自第k个分模型 \\ 0, & 否则 \end{aligned}\right. \\ j = 1,2,\cdots, N; k=1, 2,\cdots,K以上是随机变量的分布，取第1个值的概率为\alpha_1，\\取第2个值的概率为\alpha_2，……，取第K个值的概率为\alpha_K，一旦知道\gamma_1的值，就知道从第几个高斯分布中抽取y_1。\\ \begin{aligned} p(\gamma_1,y_1|\theta)= p(\gamma_1|\theta) \cdot p(y_1 | \gamma_1,\theta) \\ = \alpha^{\gamma_{11}}1 \cdot \alpha^{\gamma{12}}2 \cdots \alpha^{\gamma{1K}}K \phi(y_1|\theta_1)^{\gamma{11}} \phi(y_2|\theta_2)^{\gamma_{12}} \cdots \phi(y_1|\theta_K)^{\gamma_{1K}} = \prod_{k=1}^K [ \alpha_k \phi(y_1 | \theta_k) ]^{\gamma_{1k}} \end{aligned}\\这个是第1个样本点完全数据的密度函数。在极大化似然估计中是极大化似然函数，这需要所有样本点的联合分布，\\对于所有的样本点，概率密度函数为 \\P(y,\gamma|\theta)=\prod_{j=1}^N \prod_{k=1}^K [\alpha_k \phi(y_j | \theta_k)]^{\gamma_{jk}}\\\displaystyle \because \prod_{j=1}^N \prod_{k=1}^K \alpha_k^{\gamma_{jk}} = \prod_{k=1}^K \alpha_k^{\sum_{j=1}^N \gamma_{jk}}，\displaystyle \sum_{j=1}^N \gamma_{jk}表示在N个样本点中，一共有多少个是来自第k个高斯分布的，将该数量记为\\\displaystyle n_k=\sum_{j=1}^N \gamma_{jk}，n_1+n_2+\cdots+n_K=N\\ \displaystyle \therefore \prod_{k=1}^K \alpha_k^{\sum_{j=1}^N \gamma_{jk}} = \prod_{k=1}^K \alpha_k^{n_k}\\ \displaystyle \therefore \prod_{j=1}^N \prod_{k=1}^K \alpha_k^{\gamma_{jk}} = \prod_{k=1}^K \alpha_k^{n_k}\\ \displaystyle \therefore P(y, \gamma|\theta) = \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N \big[\phi(y_i|\theta_k)\big]^{\gamma_{jk}} = \prod_{k=1}^K \alpha_k^{n_k} \cdot \prod_{j=1}^N [ \frac{1}{\sqrt{2\pi} \sigma_k} \exp\left(-\frac{(y_j-\mu_k)^2}{2 \sigma_k^2} ) \right]^{\gamma_{jk}}\\ \displaystyle \therefore \ln P(y, \gamma|\theta) = \sum_{k=1}^K \{ n_k \ln \alpha_k + \sum_{j=1}^N \gamma_{jk} [\ln (\frac{1}{\sqrt{2\pi}}) - \ln \sigma_k - \frac{1}{2 \sigma_k^2} (y_j - \mu_k)^2] \} \\ 根据EM算法，存在一个隐变量γ，γ表示当前的y来自的高斯分布，对于第一个观测值，有γ1​=(γ11​,γ12​,⋯,γ1K​)，其中根据书中的γjk​的定义：γjk​={1,0,​第j个观测来自第k个分模型否则​j=1,2,⋯,N;k=1,2,⋯,K以上是随机变量的分布，取第1个值的概率为α1​，取第2个值的概率为α2​，……，取第K个值的概率为αK​，一旦知道γ1​的值，就知道从第几个高斯分布中抽取y1​。p(γ1​,y1​∣θ)=p(γ1​∣θ)⋅p(y1​∣γ1​,θ)=αγ11​1⋅αγ122⋯αγ1KKϕ(y1​∣θ1​)γ11ϕ(y2​∣θ2​)γ12​⋯ϕ(y1​∣θK​)γ1K​=k=1∏K​[αk​ϕ(y1​∣θk​)]γ1k​​这个是第1个样本点完全数据的密度函数。在极大化似然估计中是极大化似然函数，这需要所有样本点的联合分布，对于所有的样本点，概率密度函数为P(y,γ∣θ)=j=1∏N​k=1∏K​[αk​ϕ(yj​∣θk​)]γjk​∵j=1∏N​k=1∏K​αkγjk​​=k=1∏K​αk∑j=1N​γjk​​，j=1∑N​γjk​表示在N个样本点中，一共有多少个是来自第k个高斯分布的，将该数量记为nk​=j=1∑N​γjk​，n1​+n2​+⋯+nK​=N∴k=1∏K​αk∑j=1N​γjk​​=k=1∏K​αknk​​∴j=1∏N​k=1∏K​αkγjk​​=k=1∏K​αknk​​∴P(y,γ∣θ)=k=1∏K​αknk​​j=1∏N​[ϕ(yi​∣θk​)]γjk​=k=1∏K​αknk​​⋅j=1∏N​[2π ​σk​1​exp(−2σk2​(yj​−μk​)2​)]γjk​∴lnP(y,γ∣θ)=k=1∑K​{nk​lnαk​+j=1∑N​γjk​[ln(2π ​1​)−lnσk​−2σk2​1​(yj​−μk​)2]}

EM算法的E步，确定Q函数

$$\begin{gathered} 将隐变量都换成期望，隐变量有{\gamma }_{jk}和n_k \\ \because E(n_k)=E\left(\sum _j{\gamma }_{jk}\right)=\sum _{j}E({\gamma }_{jk})，E({\gamma }_{jk}|{\theta }^{(i)},y)=P({\gamma }_{jk}=1|{\theta }^{(i)},y)， \\ 求解期望时，是根据上一步的{\theta }^{(i)}以及观测数据所有的y_j，需要知道{\gamma }_{jk}的分布P({\gamma }_{jk}=1|{\theta }^{(i)},y)。 \\ \begin{array}{rl}\because P({\gamma }_{jk}=1|{\theta }^{(i)},y)=& \frac{P({\gamma }_{jk}=1,y_j|{\theta }^{(i)})}{P(y_j|{\theta }^{(i)})}\\ =& \frac{P({\gamma }_{jk}=1,y_j|{\theta }^{(i)})}{\sum _{k=1}^{K}P({\gamma }_{jk}=1,y_j|{\theta }^{(i)})}\\ =& \frac{P({\gamma }_{jk}=1|{\theta }^{(i)})P(y_i|{\gamma }_{jk}=1,{\theta }^{(i)})}{\sum _{k=1}^{K}P(y_j|{\gamma }_{jk}=1,{\theta }^{(i)})P({\gamma }_{jk}=1|{\theta }^{(i)})}\end{array} \\ \because {\alpha }_k=P({\gamma }_{jk}=1|\theta ),\varphi (y_i|\theta )=P(y_i|{\gamma }_{jk}=1,\theta ) \\ \therefore E({\gamma }_{jk}|y,{\theta }^{(i)})=P({\gamma }_{jk}=1|{\theta }^{(i)},y)=\frac{{\alpha }_k\varphi (y_i|{\theta }^{(i)})}{\sum _{k=1}^K{\alpha }_k\varphi (y_i|{\theta }^{(i)})}，其中{\theta }^{(i)}=({\alpha }_k^{(i)},{\theta }_k^{(i)}) \\ 将{\gamma }_{jk}关于给定y和{\theta }^{(i)}的条件下的期望记为Z_k=E({\gamma }_{jk}|y,{\theta }^{(i)})， \\ 因为各个样本之间是独立同分布的，所以Z_k是和j无关的。 \\ \therefore Q(\theta ,{\theta }^{(i)})=E_Z[lnP(y,\gamma |{\theta }^{(i)})]=\sum _{k=1}^K\left\{(N{Z}_k)\ln {\alpha }_k+Z_k\sum _{j=1}^N[\ln (\frac{1}{\sqrt{2\pi }})-\ln {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2]\right\} \end{gathered}$$

确定EM算法的M步

$$\begin{gathered} 需要估计的变量有{\alpha }_k,{\sigma }_k,{\mu }_k，然后求偏导等于0： \\ \begin{array}{l}\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\mu }_k}=0\\ \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\sigma }_k^2}=0\\ \{\begin{array}{l}\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\alpha }_k}=0\\ \sum {\alpha }_k=1\end{array}\end{array} \\ 根据上述公式可以推导出：\begin{array}{l}{\mu }_k^{(i+1)}=\frac{\sum _{j=1}^N\widehat{{\gamma }_{jk}}{y}_j}{\sum _{j=1}^N\widehat{{\gamma }_{jk}}}\\ ({\sigma }_k^2{)}^{(i+1)}=\frac{\sum _{j=1}^N\widehat{{\gamma }_{jk}}(y_i-{\mu }_k{)}^2}{\sum _{j=1}^N\widehat{{\gamma }_{jk}}}\\ {\alpha }_k^{(i+1)}=\frac{n_k}{N}=\frac{\sum _{j=1}^N\widehat{{\gamma }_{jk}}}{N}\end{array},其中\widehat{{\gamma }_{jk}}=E{\gamma }_{jk},n_k=\sum _{j=1}^{N}E{\gamma }_{jk},k=1,2,\cdots ,K \end{gathered}$$

EM算法的推广

GEM算法

$$\begin{gathered} 输入：观测数据，Q函数 \\ 输出：模型参数 \\ (1)初始化参数{\theta }^{(0)}=({\theta }_1^{(0)},{\theta }_2^{(0)},\cdots ,{\theta }_d^{(0)})，开始迭代； \\ (2)第i+1次迭代，第1步：记{\theta }^{(i)}=({\theta }_1^{(i)},{\theta }_2^{(i)},\cdots ,{\theta }_d^{(i)})为参数\theta =({\theta }_1,{\theta }_2,\cdots ,{\theta }_d)的估计值，计算 \\ \begin{array}{rlr}Q(\theta ,{\theta }^{(i)})=& E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]=& \sum _{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )\end{array} \\ (3)第2步：进行d次条件极大化：首先，在{\theta }_2^{(i)},{\theta }_3^{(i)},\cdots ,{\theta }_k^{(i)}保持不变的条件下求使得Q(\theta ,{\theta }^{(i)})达到极大的{\theta }_1^{(i+1)}； \\ 然后，在{\theta }_1={\theta }_1^{(i+1)},{\theta }_j={\theta }_j^{(j)},j=3,4,\cdots ,k的条件下求使Q(\theta ,{\theta }^{(i)})达到极大的{\theta }^{(i+1)}； \\ 如此继续，经过d次条件极大化，得到{\theta }^{(i+1)}=({\theta }_1^{(i+1)},{\theta }_2^{(i+1)},\cdots ,{\theta }_d^{(i+1)})使得Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)}) \\ (4)重复(2)和(3)，直到收敛。 \end{gathered}$$