【统计学习方法】学习笔记——EM算法及其推广


EM算法是一种迭代算法,用于含有 隐变量(hidden variable)的概率模型参数的极大似然估计,或极大后验概率估计。
EM算法的每次迭代由两步组成:E步,求期望(expectation);M步,求极大(maximization)。所以这一算法称为 期望极大算法(expectation maximization algorithm),简称EM算法

1. EM算法的引入

概率模型有时既含有观测变量( observable variable),又含有隐变量或潜在变量( (latent variable)。如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或者贝叶斯估计法估计模型参数。但是当模型含有隐变量时,就不能用这些估计方法。

EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。

1.1 EM算法

将观测数据表示为 Y = ( Y 1 , Y 2 , ⋯   , Y n ) T Y=\left(Y_{1}, Y_{2}, \cdots, Y_{n}\right)^{\mathrm{T}} Y=(Y1,Y2,,Yn)T , 未观测数据表示为 Z = ( Z 1 , Z 2 , ⋯   , Z n ) T Z=\left(Z_{1}, Z_{2}, \cdots, Z_{n}\right)^{\mathrm{T}} Z=(Z1,Z2,,Zn)T,则观测数据的似然函数为
P ( Y ∣ θ ) = ∑ Z P ( Z ∣ θ ) P ( Y ∣ Z , θ ) P(Y \mid \theta)=\sum_{Z} P(Z \mid \theta) P(Y \mid Z, \theta) P(Yθ)=ZP(Zθ)P(YZ,θ)

P ( Y ∣ θ ) = ∏ j = 1 n [ π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j ] P(Y \mid \theta)=\prod_{j=1}^{n}\left[\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}\right] P(Yθ)=j=1n[πpyj(1p)1yj+(1π)qyj(1q)1yj]
考虑求模型参数 θ = ( π , p , q ) \theta=(\pi, p, q) θ=(π,p,q)的极大似然估计, 即
θ ^ = arg ⁡ max ⁡ θ log ⁡ P ( Y ∣ θ ) \hat{\theta}=\arg \max _{\theta} \log P(Y \mid \theta) θ^=argθmaxlogP(Yθ)

这个问题没有解析解,只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。下面给出针对以上问题的 EM 算法
EM 算法首先选取参数的初值, 记作 θ ( 0 ) = ( π ( 0 ) , p ( 0 ) , q ( 0 ) ) \theta^{(0)}=\left(\pi^{(0)}, p^{(0)}, q^{(0)}\right) θ(0)=(π(0),p(0),q(0)), 然后通过下面的步骤迭代计算参数的估计值, 直至收敛为止。第 i i i次迭代参数的估计值为 θ ( i ) = ( π ( i ) , p ( i ) , q ( i ) ) \theta^{(i)}=\left(\pi^{(i)}, p^{(i)}, q^{(i)}\right) θ(i)=(π(i),p(i),q(i)) 。 EM 算法的第 i + 1 i+1 i+1 次迭代如下。
E \mathrm{E} E步: 计算在模型参数 π ( i ) , p ( i ) , q ( i ) \pi^{(i)}, p^{(i)}, q^{(i)} π(i),p(i),q(i)下观测数据 y j y_{j} yj来自掷硬币 B \mathrm{B} B 的概率
μ j ( i + 1 ) = π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j + ( 1 − π ( i ) ) ( q ( i ) ) y j ( 1 − q ( i ) ) 1 − y j \mu_{j}^{(i+1)}=\frac{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}}{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}+\left(1-\pi^{(i)}\right)\left(q^{(i)}\right)^{y_{j}}\left(1-q^{(i)}\right)^{1-y_{j}}} μj(i+1)=π(i)(p(i))yj(1p(i))1yj+(1π(i))(q(i))yj(1q(i))1yjπ(i)(p(i))yj(1p(i))1yj
M M M步: 计算模型参数的新估计值
π ( i + 1 ) = 1 n ∑ j = 1 n μ j ( i + 1 ) \pi^{(i+1)}=\frac{1}{n} \sum_{j=1}^{n} \mu_{j}^{(i+1)} π(i+1)=n1j=1nμj(i+1)
p ( i + 1 ) = ∑ j = 1 n μ j ( i + 1 ) y j ∑ j = 1 n μ 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)}} p(i+1)=j=1nμj(i+1)j=1nμj(i+1)yj
q ( i + 1 ) = ∑ j = 1 n ( 1 − μ j ( i + 1 ) y j ) ∑ j = 1 n ( 1 − μ 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)})} q(i+1)=j=1n(1μj(i+1))j=1n(1μj(i+1)yj)

EM算法与初值的选择有关,选择不同的初值可能得到不同的参数估计值。

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

假设给定观测数据 Y Y Y, 其概率分布是 P ( Y ∣ θ ) P(Y \mid \theta) P(Yθ), 其中 θ \theta θ 是需要估计的模型参数, 那么不完全数据 Y Y Y 的似然函数是 P ( Y ∣ θ ) P(Y \mid \theta) P(Yθ), 对数似然函数 L ( θ ) = log ⁡ P ( Y ∣ θ ) L(\theta)=\log P(Y \mid \theta) L(θ)=logP(Yθ) ;假设 Y Y Y Z Z Z 的联合概率分布是 P ( Y , Z ∣ θ ) P(Y, Z \mid \theta) P(Y,Zθ), 那么完全数据的对数似然函数是 log ⁡ P ( Y , Z ∣ θ ) \log P(Y, Z \mid \theta) logP(Y,Zθ)

E M \mathrm{EM} EM 算法通过迭代求 L ( θ ) = log ⁡ P ( Y ∣ θ ) L(\theta)=\log P(Y \mid \theta) L(θ)=logP(Yθ)的极大似然估计。每次迭代包含两步: E \mathrm{E} E步, 求期望; M \mathrm{M} M 步, 求极大化。
算法:EM算法
输入:观测变量数据 Y Y Y,隐变量数据 Z Z Z,联合分布 P ( Y , Z ∣ θ ) P(Y, Z|\theta) P(Y,Zθ),条件分布 P ( Z ∣ Y , θ ) P(Z|Y, \theta) P(ZY,θ)
输出:模型参数 θ \theta θ
(1)选择参数的初值 θ ( 0 ) \theta^{(0)} θ(0),开始迭代
(2)E步:记 θ ( i ) \theta^{(i)} θ(i)为第 i i i次迭代参数 θ \theta θ的估计值,在第 i + 1 i+1 i+1次迭代的 E E E步,计算 Q ( θ , θ ( i ) ) = E z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z log ⁡ P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) 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)}) Q(θ,θ(i))=Ez[logP(Y,Zθ)Y,θ(i)]=ZlogP(Y,Zθ)P(ZY,θ(i))
这里, P ( Z ∣ Y , θ ( i ) ) P(Z|Y, \theta^{(i)}) P(ZY,θ(i))是在给定观测数据 Y Y Y和当前的参数估计 θ ( i ) \theta^{(i)} θ(i)下隐变量数据 Z Z Z的条件概率分布。
(3)M步:求使 Q ( θ , θ ( i ) ) Q(\theta, \theta^{(i)}) Q(θ,θ(i))极大化的 θ \theta θ,确定第 i + 1 i+1 i+1次迭代的参数的估计值 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)
θ ( i + 1 ) = a r g max ⁡ θ Q ( θ , θ ( i ) ) \theta^{(i+1)}=arg\max_{\theta} Q(\theta, \theta^{(i)}) θ(i+1)=argθmaxQ(θ,θ(i))
(4)重复第2步和第3步,直到收敛。

函数 Q ( θ , θ ( i ) ) Q(\theta, \theta^{(i)}) Q(θ,θ(i))是EM算法的核心,称为Q函数(Q function)
定义9.1 (Q函数) 完全数据的对数似然函数 log ⁡ P ( Y , Z ∣ θ ) \log P(Y, Z|\theta) logP(Y,Zθ)关于在给定观测数据 Y Y Y和当前参数 θ ( i ) \theta^{(i)} θ(i)下对未观测数据 Z Z Z的条件概率分布 P ( Z ∣ Y , θ ( i ) ) P(Z|Y, \theta^{(i)}) P(ZY,θ(i))的期望称为Q函数,即 Q ( θ , θ ( i ) ) = E z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] Q(\theta, \theta^{(i)})=E_z[\log P(Y, Z|\theta)|Y, \theta^{(i)}] Q(θ,θ(i))=Ez[logP(Y,Zθ)Y,θ(i)]

下面关于 EM 算法作几点说明:

  • 步骤 (1)参数的初值可以任意选择, 但需注意 E M \mathrm{EM} EM算法对初值是敏感的。
  • 步骤 (2) E \mathrm{E} E步求 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i)) 。Q 函数式中 Z 是未观测数据, Y 是观测数据。注 意 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))的第 1 个变元表示要极大化的参数, 第 2 个变元表示参数的当前估计值。每次迭代实际在求 Q 函数及其极大。
  • 步骤 (3) M \mathrm{M} M步求 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))的极大化, 得到 θ ( i + 1 ) \theta^{(i+1)} θ(i+1), 完成一次迭代 θ ( i ) → θ ( i + 1 ) \theta^{(i)} \rightarrow \theta^{(i+1)} θ(i)θ(i+1)。 后面将证明每次迭代使似然函数增大或达到局部极值。
  • 步骤 (4)给出停止迭代的条件, 一般是对较小的正数 ε 1 , ε 2 \varepsilon_{1}, \varepsilon_{2} ε1,ε2, 若满足
    ∥ θ ( i + 1 ) − θ ( i ) ∥ < ε 1  或  ∥ Q ( θ ( i + 1 ) , θ ( i ) ) − Q ( θ ( i ) , θ ( i ) ) ∥ < ε 2 \left\|\theta^{(i+1)}-\theta^{(i)}\right\|<\varepsilon_{1} \quad \text { 或 }\left\|Q\left(\theta^{(i+1)}, \theta^{(i)}\right)-Q\left(\theta^{(i)}, \theta^{(i)}\right)\right\|<\varepsilon_{2} θ(i+1)θ(i)<ε1  Q(θ(i+1),θ(i))Q(θ(i),θ(i))<ε2
    则停止迭代。备注:实际工程中可能会设置epoch迭代停止

1.2 EM算法的导出

EM算法本质是估计一个密度函数,估计密度函数通常情况下是通过观测值,最大化似然函数去估计参数,而EM算法针对含有隐变量的密度函数估计,此时直接最大化似然函数会比较困难。

待估计参数 θ \theta θ, 观测数据: Y = { y 1 , y 2 , . . . , y N } Y=\{y_1, y_2,...,y_N\} Y={y1,y2,...,yN} ,隐变量: Z = { z 1 , z 2 , . . . , N } Z=\{z_1, z_2,...,_N\} Z={z1,z2,...,N} ,似然函数: L ( θ ) = log ⁡ P ( Y ∣ θ ) = log ⁡ ∑ Z P ( Y , Z ∣ θ ) = log ⁡ ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) L(\theta) = \log P(Y|\theta) = \log \sum_Z P(Y,Z|\theta) = \log (\sum_Z P(Y|Z,\theta)P(Z|\theta)) L(θ)=logP(Yθ)=logZP(Y,Zθ)=log(ZP(YZ,θ)P(Zθ))

对于给定 θ \theta θ一个初值,然后进行迭代逐渐最大化似然函数,
L ( θ ( i + 1 ) ) − L ( θ ( i ) ) = log ⁡ ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) − log ⁡ P ( Y ∣ θ ( i ) ) L(\theta^{(i+1)})-L(\theta^{(i)})=\log (\sum_Z P(Y|Z,\theta)P(Z|\theta))-\log P(Y|\theta^{(i)}) L(θ(i+1))L(θ(i))=log(ZP(YZ,θ)P(Zθ))logP(Yθ(i))
由Jensen不等式: f ( ∑ i α i x i ) ≥ ∑ i α i f ( x i ) ⇒ log ⁡ E ( Y ) ≥ E ( log ⁡ ( Y ) ) f(\sum_i \alpha_i x_i)\ge \sum_i \alpha_i f(x_i) \Rightarrow \log E(Y) \ge E(\log (Y)) f(iαixi)iαif(xi)logE(Y)E(log(Y)),其中 0 ≤ α i ≤ 1 0\le \alpha_i \le 1 0αi1
L ( θ ( i + 1 ) ) − L ( θ ( i ) ) = log ⁡ ∑ Z P ( Y ∣ Z , θ ( i ) ) P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ( i ) ) − log ⁡ P ( Y ∣ θ ( i ) ) ≥ ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) − log ⁡ P ( Y ∣ θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) L(\theta^{(i+1)})-L(\theta^{(i)}) = \log{\sum_Z P(Y|Z,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta^{(i)})}}-\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)})P(Y|\theta^{(i)})}} L(θ(i+1))L(θ(i))=logZP(YZ,θ(i))P(YZ,θ(i))P(YZ,θ)P(Zθ)logP(Yθ(i))ZP(ZY,θ(i))logP(ZY,θ(i))P(YZ,θ)P(Zθ)logP(Yθ(i))=ZP(ZY,θ(i))logP(ZY,θ(i))P(Yθ(i))P(YZ,θ)P(Zθ)
B ( θ , θ ( i ) ) = L ( θ ( i ) ) + ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) B(\theta, \theta^{(i)})=L(\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)})}} B(θ,θ(i))=L(θ(i))+ZP(ZY,θ(i))logP(ZY,θ(i))P(Yθ(i))P(YZ,θ)P(Zθ)
则有: L ( θ ) ≥ B ( θ , θ ( i ) ) L(\theta) \ge B(\theta, \theta^{(i)}) L(θ)B(θ,θ(i))
即函数 B ( θ , θ ( i ) ) B(\theta, \theta^{(i)}) B(θ,θ(i)) L ( θ ) L(\theta) L(θ)的一个下界,而且由上式可知 L ( θ ( i ) ) = B ( θ ( i ) , θ ( i ) ) L(\theta^{(i)}) = B(\theta^{(i)}, \theta^{(i)}) L(θ(i))=B(θ(i),θ(i))
因此,任何可以使 B ( θ ( i ) , θ ( i ) ) B(\theta^{(i)}, \theta^{(i)}) B(θ(i),θ(i))增大的 θ \theta θ,也可以使 L ( θ ) L(\theta) L(θ)增大。选择 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)使得 B ( θ , θ ( i ) ) B(\theta, \theta^{(i)}) B(θ,θ(i))达到极大,即
θ ( i + 1 ) = a r g max ⁡ θ B ( θ ( i ) , θ ( i ) ) \theta^{(i+1)}=arg \max_{\theta} B(\theta^{(i)}, \theta^{(i)}) θ(i+1)=argθmaxB(θ(i),θ(i))
去除跟 θ \theta θ无关的项,有:
θ ( i + 1 ) = a r g max ⁡ θ ( ∑ Z P ( Z ∣ Y , θ ( i ) ) ( log ⁡ P ( Y ∣ Z , θ ) ) P ( Z ∣ θ ) ) = a r g max ⁡ θ ∑ Z P ( Z ∣ Y , θ ( i ) ) ( log ⁡ P ( Y , Z ∣ θ ) ) \theta^{(i+1)}=arg \max_{\theta} (\sum_Z P(Z|Y,\theta^{(i)}) (\log P(Y|Z, \theta))P(Z|\theta))\\\\ = arg\max_\theta \sum_Z P(Z|Y, \theta^{(i)}) (\log P(Y, Z|\theta)) θ(i+1)=argθmax(ZP(ZY,θ(i))(logP(YZ,θ))P(Zθ))=argθmaxZP(ZY,θ(i))(logP(Y,Zθ))
Q ( θ , θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) ( log ⁡ P ( Y , Z ∣ θ ) ) Q\left(\theta, \theta^{(i)}\right)=\sum_Z P(Z|Y, \theta^{(i)}) (\log P(Y, Z|\theta)) Q(θ,θ(i))=ZP(ZY,θ(i))(logP(Y,Zθ))
Q i ( Z ( i ) ) : = P ( Z ∣ Y , θ ( i ) ) Q_i(Z^{(i)}):=P(Z|Y, \theta^{(i)}) Qi(Z(i)):=P(ZY,θ(i))就是EM算法中Expectation,
a r g max ⁡ θ Q ( θ , θ ( i ) ) arg \max_\theta Q\left(\theta, \theta^{(i)}\right) argmaxθQ(θ,θ(i))就是EM算法中Maximization。

EM 算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
EM算法的直观理解

1.3 EM算法在非监督学习中的应用

监督学习是由训练数据 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x N , y N ) } \left\{\left(x_{1}, y_{1}\right),\left(x_{2}, y_{2}\right), \cdots,\left(x_{N}, y_{N}\right)\right\} {(x1,y1),(x2,y2),,(xN,yN)}学习条件概率分布 P ( Y ∣ X ) P(Y \mid X) P(YX)或决策函数 Y = f ( X ) Y=f(X) Y=f(X)作为模型, 用于分类、回归、标注等任务。这时训练数据中的每个样本点由输入和输出对组成。
有时训练数据只有输入没有对应的输出 { ( x 1 , ∙ ) , ( x 2 , ∙ ) , ⋯   , ( x N , ∙ ) } \left\{\left(x_{1}, \bullet\right),\left(x_{2}, \bullet\right), \cdots,\left(x_{N}, \bullet\right)\right\} {(x1,),(x2,),,(xN,)}, 从这样的数据学习模型称为无监督学习问题

EM 算法可以用于生成模型的无监督学习。生成模型由联合概率分布 P ( X , Y ) P(X,Y) P(X,Y) 表示, 可以认为无监督学习训练数据是联合概率分布产生的数据。 X X X 为观测数据, Y Y Y 为未观测数据。

2. EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。EM最大的优点就是简单性和普适性。
定理9.1 设 P ( Y ∣ θ ) P(Y \mid \theta) P(Yθ) 为观测数据的似然函数, θ ( i ) ( i = 1 , 2 , ⋯   ) \theta^{(i)}(i=1,2, \cdots) θ(i)(i=1,2,) E M \mathrm{EM} EM 算法得到的参数估计序列, P ( Y ∣ θ ( i ) ) ( i = 1 , 2 , ⋯   ) P\left(Y \mid \theta^{(i)}\right)(i=1,2, \cdots) P(Yθ(i))(i=1,2,)为对应的似然函数序列, 则 P ( Y ∣ θ ( i ) ) P\left(Y \mid \theta^{(i)}\right) P(Yθ(i))是单调递增的,即
P ( Y ∣ θ ( i + 1 ) ) ⩾ P ( Y ∣ θ ( i ) ) P\left(Y \mid \theta^{(i+1)}\right) \geqslant P\left(Y \mid \theta^{(i)}\right) P(Yθ(i+1))P(Yθ(i))

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

在许多情况下,EM算法是学习高斯混合模型(Gaussian misture model)

3.1 高斯混合模型

定义9.2 (高斯混合模型) 高斯混合模型是指具有如下形式的概率分布模型:
P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y|\theta) = \sum_{k=1}^K \alpha_k \phi(y|\theta_k) P(yθ)=k=1Kαkϕ(yθk)
其中, α k \alpha_k αk是系数, α k ≥ 0 , ∑ k = 1 K α k = 1 ; ϕ ( y ∣ θ k ) \alpha_k \ge 0, \sum_{k=1}^K \alpha_k=1; \phi(y|\theta_k) αk0,k=1Kαk=1;ϕ(yθk)是高斯分布密度, θ k = ( μ k , θ k 2 ) \theta_k=(\mu_k, \theta_k^2) θk=(μkθk2)
ϕ ( y ∣ θ k ) = 1 2 π σ k e x p ( − ( y − μ k ) 2 2 σ k 2 ) \phi (y|\theta_k)=\frac{1}{\sqrt{2\pi} \sigma_k}exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2}) ϕ(yθk)=2π σk1exp(2σk2(yμk)2)
称为第 k k k个分模型。
一般混合模型可以由任意概率分布密度代替上式中的高斯分布密度。

3.2 高斯混合模型参数估计的EM算法

观测数据 y 1 , y 2 , ⋯   , y N y_{1}, y_{2}, \cdots, y_{N} y1,y2,,yN由高斯混合模型生成,
P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y \mid \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y \mid \theta_{k}\right) P(yθ)=k=1Kαkϕ(yθk)
其中, θ = ( α 1 , α 2 , ⋯   , α K ; θ 1 , θ 2 , ⋯   , θ K ) \theta=\left(\alpha_{1}, \alpha_{2}, \cdots, \alpha_{K} ; \theta_{1}, \theta_{2}, \cdots, \theta_{K}\right) θ=(α1,α2,,αK;θ1,θ2,,θK) 。用 EM算法估计高斯混合模型的参数 θ \theta θ

  1. 明确隐变量,写出完全数据的对数似然函数
    可以设想观测数据 y j , j = 1 , 2 , ⋯   , N y_{j}, j=1,2, \cdots, N yj,j=1,2,,N, 是这样产生的: 首先依概率 α k \alpha_{k} αk选择第 k k k个高斯分布分模型 ϕ ( y ∣ θ k ) \phi\left(y \mid \theta_{k}\right) ϕ(yθk), 然后依第 k k k 个分模型的概率分布 ϕ ( y ∣ θ k ) \phi\left(y \mid \theta_{k}\right) ϕ(yθk)生成观测数据 y j y_{j} yj 。 这时观测数据 y j , j = 1 , 2 , ⋯   , N y_{j}, j=1,2, \cdots, N yj,j=1,2,,N, 是已知的; 反映观测数据 y j y_{j} yj来自第 k k k个分模型的数据是未知的, k = 1 , 2 , ⋯   , K k=1,2, \cdots, K k=1,2,,K, 以隐变量 γ j k \gamma_{j k} γjk表示, 其定义如下:
    γ j k = { 1 第 j 个 观 测 来 自 第 k 个 分 模 型 0 否 则 \gamma_{jk}=\left\{\begin{matrix} 1 & 第j个观测来自第k个分模型\\ 0 & 否则 \end{matrix}\right. γjk={10jk
    j = 1 , 2 , . . . , N ; k = 1 , 2 , . . . , K j=1,2,..., N; k=1,2,...,K j=1,2,...,N;k=1,2,...,K
    γ j k \gamma_{j k} γjk 0 − 1 0-1 01随机变量。
    有了观测数据 y j y_{j} yj及未观测数据 γ j k \gamma_{j k} γjk, 那么完全数据是
    ( y j , γ j 1 , γ j 2 , ⋯   , γ j K ) , j = 1 , 2 , ⋯   , N \left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}\right), \quad j=1,2, \cdots, N (yj,γj1,γj2,,γjK),j=1,2,,N
    完全数据的似然函数:
    P ( y , γ ∣ θ ) = ∏ j = 1 N P ( y j , γ j 1 , γ j 2 , ⋯   , γ j K ∣ θ ) = ∏ k = 1 K ∏ j = 1 N [ α k ϕ ( y j ∣ θ k ) ] γ j k = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y j ∣ θ k ) ] γ j k ∏ k = 1 K α k n k ∏ j = 1 N [ 1 2 π σ k e x p ( − ( y j − μ k ) 2 2 σ k 2 ) ] γ j k P(y,\gamma |\theta)=\prod_{j=1}^N P(y_j, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j K}|\theta)\\\\ = \prod_{k=1}^K \prod_{j=1}^N [\alpha_k \phi(y_j|\theta_k)]^{\gamma_{j k}}\\\\ =\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\phi(y_j|\theta_k)]^{\gamma_{j k}}\\\\ \prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\frac{1}{\sqrt{2\pi} \sigma_k}exp(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2})]^{\gamma_{j k}} P(y,γθ)=j=1NP(yj,γj1,γj2,,γjKθ)=k=1Kj=1N[αkϕ(yjθk)]γjk=k=1Kαknkj=1N[ϕ(yjθk)]γjkk=1Kαknkj=1N[2π σk1exp(2σk2(yjμk)2)]γjk
    式中, n k = ∑ j = 1 N γ j k , ∑ k = 1 K n k = N n_k=\sum_{j=1}^N\gamma_{j k}, \sum_{k=1}^K n_k=N nk=j=1Nγjk,k=1Knk=N
    那么, 完全数据的对数似然函数为
    log ⁡ P ( y , γ ∣ θ ) = ∑ k = 1 K { n k log ⁡ α k + ∑ j = 1 N γ j k [ log ⁡ ( 1 2 π ) − log ⁡ σ k − 1 2 σ k 2 ( y j − μ k ) 2 ] } \log P(y, \gamma \mid \theta)=\sum_{k=1}^{K}\left\{n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\} logP(y,γθ)=k=1K{nklogαk+j=1Nγjk[log(2π 1)logσk2σk21(yjμk)2]}

  2. EM算法的E步:确定Q函数
    确定Q函数

  3. 确定EM算法的M步
    迭代的 M \mathrm{M} M步是求函数 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i)) θ \theta θ的极大值, 即求新一轮迭代的模型参数:
    θ ( i + 1 ) = arg ⁡ max ⁡ θ Q ( θ , θ ( i ) ) \theta^{(i+1)}=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right) θ(i+1)=argθmaxQ(θ,θ(i))
    μ ^ k \hat{\mu}_{k} μ^k, σ ^ k 2 \hat{\sigma}_{k}^{2} σ^k2 α ^ k , k = 1 , 2 , ⋯   , K \hat{\alpha}_{k}, k=1,2, \cdots, K α^k,k=1,2,,K 表示 θ ( i + 1 ) \theta^{(i+1)} θ(i+1) 的各参数。求
    μ ^ k , σ ^ k 2 \hat{\mu}_{k}, \hat{\sigma}_{k}^{2} μ^k,σ^k2 只需将式④分别对 μ k , σ k 2 \mu_{k}, \sigma_{k}^{2} μk,σk2 求偏导数并令其为 0, 即可得到; 求 α ^ k \hat{\alpha}_{k} α^k 是在 ∑ k = 1 K α k = 1 \sum_{k=1}^{K} \alpha_{k}=1 k=1Kαk=1条件下求偏导数并令其为 0 得到的。结果如下:
    μ k ^ = ∑ j = 1 N γ ^ j k y j ∑ j = 1 N γ ^ j k , k = 1 , 2 , . . . , K \hat{\mu_k}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}}, k=1,2,...,K μk^=j=1Nγ^jkj=1Nγ^jkyj,k=1,2,...,K
    σ ^ k 2 = ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 ∑ j = 1 N γ ^ j k , k = 1 , 2 , . . . , K \hat{\sigma}_k^2=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}}, k=1,2,...,K σ^k2=j=1Nγ^jkj=1Nγ^jk(yjμk)2,k=1,2,...,K
    α ^ k = n k N = ∑ j = 1 N γ ^ j k N , k = 1 , 2 , . . . , K \hat{\alpha}_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N}, k=1,2,...,K α^k=Nnk=Nj=1Nγ^jk,k=1,2,...,K
    重复以上计算, 直到对数似然函数值不再有明显的变化为止。

算法:高斯混合模型参数估计的EM算法
输入: 观测数据 y 1 , y 2 , ⋯   , y N y_{1}, y_{2}, \cdots, y_{N} y1,y2,,yN, 高斯混合模型;
输出:高斯混合模型参数。
(1)取参数的初始值开始迭代
(2)E步:依据当前模型参数,计算分模型 k k k对观测数据 y j y_j yj的响应度
γ ^ j k = α k ϕ ( y j ∣ θ k ) ∑ k = 1 K α k ϕ ( y j ∣ θ k ) , j = 1 , 2 , . . . , N ; k = 1 , 2 , . . . , K \hat{\gamma}_{jk}=\frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^K \alpha_k \phi(y_j|\theta_k)},j=1,2,...,N; k=1,2,...,K γ^jk=k=1Kαkϕ(yjθk)αkϕ(yjθk)j=1,2,...,N;k=1,2,...,K
(3)M步:计算新一轮迭代的模型参数
μ k ^ = ∑ j = 1 N γ ^ j k y j ∑ j = 1 N γ ^ j k , k = 1 , 2 , . . . , K \hat{\mu_k}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}}, k=1,2,...,K μk^=j=1Nγ^jkj=1Nγ^jkyj,k=1,2,...,K
σ ^ k 2 = ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 ∑ j = 1 N γ ^ j k , k = 1 , 2 , . . . , K \hat{\sigma}_k^2=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}}, k=1,2,...,K σ^k2=j=1Nγ^jkj=1Nγ^jk(yjμk)2,k=1,2,...,K
α ^ k = ∑ j = 1 N γ ^ j k N , k = 1 , 2 , . . . , K \hat{\alpha}_k=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N}, k=1,2,...,K α^k=Nj=1Nγ^jk,k=1,2,...,K
(4)重复第(2)步和第(3)步,直到收敛。

4. EM算法的推广

EM算法还可以解释为F函数(F function)的极大—极大算法(maximization-maximization algorithm),基于这个解释有若干变形与推广,如广义期望极大(generalized expectation maximization, GEM)算法。

4.1 F函数的极大—极大算法

F 函数 假设隐变量数据 Z Z Z 的概率分布为 P ~ ( Z ) \tilde{P}(Z) P~(Z) , 定义分布 P ~ \tilde{P} P~ 与参数 θ \theta θ的函数 F ( P ~ , θ ) F(\tilde{P}, \theta) F(P~,θ)如下:
F ( P ~ , θ ) = E P ~ [ log ⁡ P ( Y , Z ∣ θ ) ] + H ( P ~ ) = F ( P ~ , θ ) = ∑ Z [ P ~ ( Z ) log ⁡ P ( Y , Z ∣ θ ) ] − ∑ Z [ P ~ ( Z ) log ⁡ P ~ ( Z ) ] F(\tilde{P}, \theta)=E_{\tilde{P}}[\log P(Y, Z \mid \theta)]+H(\tilde{P})=F(\tilde{P}, \theta)=\sum_Z[ \tilde{P}(Z)\log P(Y, Z \mid \theta)]-\sum_Z[ \tilde{P}(Z)\log \tilde{P}(Z)] F(P~,θ)=EP~[logP(Y,Zθ)]+H(P~)=F(P~,θ)=Z[P~(Z)logP(Y,Zθ)]Z[P~(Z)logP~(Z)]
称为 F函数。式中 H ( P ~ ) = − E P ~ log ⁡ P ~ ( Z ) H(\tilde{P})=-E_{\tilde{P}} \log \tilde{P}(Z) H(P~)=EP~logP~(Z)是分布 P ~ ( Z ) \tilde{P}(Z) P~(Z) 的熵。

引理1 对于固定的 θ \theta θ, 存在唯一的分布 P ~ θ \tilde{P}_{\theta} P~θ极大化 F ( P ~ , θ ) F(\tilde{P}, \theta) F(P~,θ), 这时 P ~ θ \tilde{P}_{\theta} P~θ
由下式给出:
P ~ θ ( Z ) = P ( Z ∣ Y , θ ) \tilde{P}_{\theta}(Z)=P(Z \mid Y, \theta) P~θ(Z)=P(ZY,θ)
并且 P ~ θ \tilde{P}_{\theta} P~θ θ \theta θ 连续变化。
引理2 P ~ θ ( Z ) = P ( Z ∣ Y , θ ) \tilde{P}_{\theta}(Z)=P(Z \mid Y, \theta) P~θ(Z)=P(ZY,θ) , 则 F ( P ~ , θ ) = l o g P ( Y ∣ θ ) F(\tilde{P},\theta)=logP(Y\mid\theta) F(P~,θ)=logP(Yθ)
定理3 L ( θ ) = log ⁡ P ( Y ∣ θ ) L(\theta)=\log P(Y \mid \theta) L(θ)=logP(Yθ) 为观测数据的对数似然函数, θ ( i ) , i = 1 , 2 , ⋯ \theta^{(i)}, i=1,2, \cdots θ(i),i=1,2, E M \mathrm{EM} EM 算法得到的参数估计序列, 如果 F ( P ~ , θ ) F(\tilde{P}, \theta) F(P~,θ) P ~ ∗ \tilde{P}^{*} P~ θ ∗ \theta^{*} θ有局部极大值, 那么 L ( θ ) L(\theta) L(θ) 也在 θ ∗ \theta^{*} θ有局部极大值。类似地, 如果 F ( P ~ , θ ) F(\tilde{P}, \theta) F(P~,θ) P ~ ∗ \tilde{P}^{*} P~
θ ∗ \theta^{*} θ达到全局最大值, 那么 L ( θ ) L(\theta) L(θ)也在 θ ∗ \theta^{*} θ达到全局最大值。
定理4 EM 算法的一次迭代可由 F 函数的极大-极大算法实现. 设 θ ( i ) \theta^{(i)} θ(i)为第 i i i 次迭代参数 θ \theta θ 的估计, P ~ ( i ) \tilde{P}^{(i)} P~(i) 为第 i i i 次迭代函数 P ~ \tilde{P} P~的估计。在第 i + 1 i+1 i+1次迭代的两步为:
(1) 对固定的 θ ( i ) \theta^{(i)} θ(i) P ~ ( i + 1 ) \tilde{P}^{(i+1)} P~(i+1) 使 F ( P ~ , θ ( i ) ) F\left(\tilde{P}, \theta^{(i)}\right) F(P~,θ(i))极大化;
(2) 对固定的 P ~ ( i + 1 ) \tilde{P}^{(i+1)} P~(i+1) θ ( i + 1 ) \theta^{(i+1)} θ(i+1)使 F ( P ~ ( i + 1 ) , θ ) F\left(\tilde{P}^{(i+1)}, \theta\right) F(P~(i+1),θ)极大化.

4.2 GEM算法

GEM 算法 1
输入: 观测数据, F 函数;
输出:模型参数。
(1) 初始化参数 θ ( 0 ) \theta^{(0)} θ(0),开始迭代:
(2) 第 i + 1 i+1 i+1 次迭代, 第 1 步: 记 θ ( i ) \theta^{(i)} θ(i)为参数 θ \theta θ的估计值, P ~ ( i ) \tilde{P}^{(i)} P~(i)为函数 P ~ \tilde{P} P~的估计,
P ~ ( i + 1 ) \tilde{P}^{(i+1)} P~(i+1) 使 P ~ \tilde{P} P~ 极大化 F ( P ~ , θ ( i ) ) F\left(\tilde{P}, \theta^{(i)}\right) F(P~,θ(i));
(3) 第 2 步: 求 θ ( i + 1 ) \theta^{(i+1)} θ(i+1) 使 F ( P ~ ( i + 1 ) , θ ) F\left(\tilde{P}^{(i+1)}, \theta\right) F(P~(i+1),θ)极大化;
(4) 重复
(2) 和 ( 3 ) (3)(3), 直到收敛。
在 GEM算法 1 中, 有时求 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))的极大化是很困难的。下面介绍的GEM算法 2 和GEM算法 3 并不是直接求 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)使 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))达到极大的 θ \theta θ, 而是找一个 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)使得 Q ( θ ( i + 1 ) , θ ( i ) ) > Q ( θ ( i ) , θ ( i ) ) Q\left(\theta^{(i+1)}, \theta^{(i)}\right)>Q\left(\theta^{(i)}, \theta^{(i)}\right) Q(θ(i+1),θ(i))>Q(θ(i),θ(i)) .

GEM 算法 2
输入: 观测数据, Q 函数;
输出: 模型参数。
(1) 初始化参数 θ ( 0 ) \theta^{(0)} θ(0), 开始迭代:
(2) 第 i + 1 i+1 i+1次迭代, 第 1 步: 记 θ ( i ) \theta^{(i)} θ(i)为参数 θ \theta θ 的估计值, 计算
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z P ( Z ∣ Y , θ ( i ) ) log ⁡ P ( Y , Z ∣ θ ) Q(θ,θ^{(i)})=E_Z[\log P(Y,Z∣θ)∣Y,θ^{(i)}]=\sum_ZP(Z∣Y,θ^{(i)})\log P(Y,Z∣θ) Q(θ,θ(i))=EZ[logP(Y,Zθ)Y,θ(i)]=ZP(ZY,θ(i))logP(Y,Zθ)

(3) 第 2 步: 求 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)使
Q ( θ ( i + 1 ) , θ ( i ) ) > Q ( θ ( i ) , θ ( i ) ) Q\left(\theta^{(i+1)}, \theta^{(i)}\right)>Q\left(\theta^{(i)}, \theta^{(i)}\right) Q(θ(i+1),θ(i))>Q(θ(i),θ(i))
(2) 和 ( 3 ) , 直到收敛。
(4) 重复
当参数 θ \theta θ的维数为 d ( d ⩾ 2 ) d(d \geqslant 2) d(d2) 时, 可采用一种特殊的GEM算法, 它将EM算法的M 步分解为 d 次条件极大化, 每次只改变参数向量的一个分量, 其余分量不改变。

GEM 算法 3
输入: 观测数据, Q 函数;
输出: 模型参数。
(1) 初始化参数 θ ( 0 ) = ( θ 1 ( 0 ) , θ 2 ( 0 ) , ⋯   , θ d ( 0 ) ) \theta^{(0)}=\left(\theta_{1}^{(0)}, \theta_{2}^{(0)}, \cdots, \theta_{d}^{(0)}\right) θ(0)=(θ1(0),θ2(0),,θd(0)), 开始迭代;
(2) 第 i + 1 i+1 i+1 次迭代, 第 1 步: 记 θ ( i ) = ( θ 1 ( i ) , θ 2 ( i ) , ⋯   , θ d ( i ) ) \theta^{(i)}=\left(\theta_{1}^{(i)}, \theta_{2}^{(i)}, \cdots, \theta_{d}^{(i)}\right) θ(i)=(θ1(i),θ2(i),,θd(i))为参数 θ = ( θ 1 , θ 2 , ⋯   , θ d ) \theta=\left(\theta_{1}, \theta_{2}, \cdots, \theta_{d}\right) θ=(θ1,θ2,,θd)的估计值, 计算
Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = ∑ Z P ( Z ∣ y , θ ( i ) ) log ⁡ P ( Y , Z ∣ θ ) Q(θ,θ^{(i)})=E_Z[\log P(Y,Z∣θ)∣Y,θ^{(i)}]=\sum_Z P(Z∣y,θ^{(i)})\log P(Y,Z∣θ) Q(θ,θ(i))=EZ[logP(Y,Zθ)Y,θ(i)]=ZP(Zy,θ(i))logP(Y,Zθ)
(3) 第 2 步: 进行 d 次条件极大化:
首先, 在 θ 2 ( i ) , ⋯   , θ d ( i ) \theta_{2}^{(i)}, \cdots, \theta_{d}^{(i)} θ2(i),,θd(i)保持不变的条件下求使 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))达到极大的 θ 1 ( i + 1 ) \theta_{1}^{(i+1)} θ1(i+1);
然后, 在 θ 1 = θ 1 ( i + 1 ) , θ j = θ j ( i ) , j = 3 , 4 , ⋯   , d \theta_{1}=\theta_{1}^{(i+1)}, \theta_{j}=\theta_{j}^{(i)}, j=3,4, \cdots, d θ1=θ1(i+1),θj=θj(i),j=3,4,,d的条件下求使 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))达到极大的 θ 2 ( i + 1 ) \theta_{2}^{(i+1)} θ2(i+1);
如此继续, 经过 d 次条件极大化, 得到 θ ( i + 1 ) = ( θ 1 ( i + 1 ) , θ 2 ( i + 1 ) , ⋯   , θ d ( i + 1 ) ) \theta^{(i+1)}=\left(\theta_{1}^{(i+1)}, \theta_{2}^{(i+1)}, \cdots, \theta_{d}^{(i+1)}\right) θ(i+1)=(θ1(i+1),θ2(i+1),,θd(i+1))使得
Q ( θ ( i + 1 ) , θ ( i ) ) > Q ( θ ( i ) , θ ( i ) ) Q\left(\theta^{(i+1)}, \theta^{(i)}\right)>Q\left(\theta^{(i)}, \theta^{(i)}\right) Q(θ(i+1),θ(i))>Q(θ(i),θ(i))
(4)重复(2) 和 ( 3 ), 直到收敛。

总结

  • EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。
  • 在构建具体的EM算法时,重要的是定义Q函数,每次迭代中,EM算法通过极大化Q函数来增大对数似然函数 L ( θ ) L(\theta) L(θ)
  • EM算法的结果和初值的选取有很大关系,不同的初值会得到不同的解,有可能收敛到局部最优。

内容来源

[1] 《统计学习方法》 李航
[2] https://blog.csdn.net/wl1780852311/article/details/119613563
[3] https://zhuanlan.zhihu.com/p/126320789

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

镰刀韭菜

看在我不断努力的份上,支持我吧

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值