EM算法
极大似然估计
极大似然估计:(maximum likelihood estimate, MLE) 是一种常用的模型参数估计方法。它假设观测样本出现的概率最大,也即样本联合概率(也称似然函数)取得最大值。
为求解方便,对样本联合概率取对数似然函数
log L ( θ ) = log P ( X ∣ θ ) = ∑ i = 1 N log P ( x i ∣ θ ) \log L(\theta) =\log\mathbb P(X|\theta)=\sum_{i=1}^N\log \mathbb P(\mathbf x_i|\theta) logL(θ)=logP(X∣θ)=i=1∑NlogP(xi∣θ)
优化目标是最大化对数似然函数
θ ^ = arg max θ ∑ i = 1 N log P ( x i ∣ θ ) \hat\theta=\arg\max_{\theta}\sum_{i=1}^N\log \mathbb P(\mathbf x_i|\theta) θ^=argθmaxi=1∑NlogP(xi∣θ)
假设瓜田里有两种类型的西瓜🍉,瓜农随机抽取了10个西瓜,来了解西瓜的重量分布 p ( x ∣ θ ) p(x|\theta) p(x∣θ),记录结果如下:
变量 | 样本 |
---|---|
西瓜重量 x x x | 5.3 , 5.7, 4.7, 4.3, 3.2, 4.9, 4.1, 3.5, 3.8, 1.7 |
西瓜品种 z z z | 1, 1, 1, 1, 2, 2, 2, 2, 2, 2 |
其中,西瓜的品种 z z z 是离散分布 P ( z = k ) = π k \mathbb P(z=k)=\pi_k P(z=k)=πk,一般假设两种类型的西瓜服从均值和方差不同的高斯分布 N ( μ 1 , σ 1 2 ) N(\mu_1,\sigma^2_1) N(μ1,σ12)和 N ( μ 2 , σ 2 2 ) N(\mu_2,\sigma^2_2) N(μ2,σ22)。由全概率公式,西瓜重量的概率密度模型
p ( x ; θ ) = π 1 N ( x ; μ 1 , σ 1 2 ) + π 2 N ( x ; μ 2 , σ 2 2 ) p(x;\theta)=\pi_1\mathcal N(x;\mu_1,\sigma^2_1)+\pi_2\mathcal N(x;\mu_2,\sigma^2_2) p(x;θ)=π1N(x;μ1,σ12)+π2N(x;μ2,σ22)
我们尝试用极大似然估计求解参数 θ = ( π 1 , π 2 , μ 1 , σ 1 2 , μ 2 , σ 2 2 ) \theta=(\pi_1,\pi_2,\mu_1,\sigma^2_1,\mu_2,\sigma^2_2) θ=(π1,π2,μ1,σ12,μ2,σ22)。
优化目标函数
max θ ∑ z i = 1 log π 1 N ( x i ; μ 1 , σ 1 2 ) + ∑ z i = 2 log π 2 N ( x i ; μ 2 , σ 2 2 ) s.t. π 1 + π 2 = 1 \max_{\theta}\sum_{z_i=1}\log \pi_1\mathcal N(x_i;\mu_1,\sigma_1^2)+\sum_{z_i=2}\log \pi_2\mathcal N(x_i;\mu_2,\sigma_2^2) \\ \text{s.t. } \pi_1+\pi_2=1 θmaxzi=1∑logπ1N(xi;μ1,σ12)+zi=2∑logπ2N(xi;μ2,σ22)s.t. π1+π2=1
使用拉格朗日乘子法容易求得
π 1 = 0.4 , π 2 = 0.6 μ 1 = 5 , σ 1 2 = 0.5 4 2 μ 2 = 3.53 , σ 2 2 = 0.9 8 2 \pi_1=0.4,\quad \pi_2=0.6 \\ \mu_1=5,\quad \sigma_1^2=0.54^2 \\ \mu_2=3.53,\quad \sigma_2^2=0.98^2 \\ π1=0.4,π2=0.6μ1=5,σ12=0.542μ2=3.53,σ22=0.982
最终得到
p ( x ) = 0.4 × N ( x ; 5 , 0.5 4 2 ) + 0.6 × N ( x ; 3.53 , 0.9 8 2 ) p(x)=0.4\times\mathcal N(x;5,0.54^2)+0.6\times\mathcal N(x;3.53,0.98^2) p(x)=0.4×N(x;5,0.542)+0.6×N(x;3.53,0.982)
但是,实际中如果瓜农无法辩识标记西瓜的品种,此时概率分布函数变为
p ( x ; θ ) = π N ( x ; μ 1 , σ 1 2 ) + ( 1 − π ) N ( x ; μ 2 , σ 2 2 ) p(x;\theta)=\pi\mathcal N(x;\mu_1,\sigma^2_1)+(1-\pi)\mathcal N(x;\mu_2,\sigma^2_2) p(x;θ)=πN(x;μ1,σ12)+(1−π)N(x;μ2,σ22)
其中品种 z z z 成为隐藏变量。对数似然函数变为
log L ( θ ) = ∑ i log ( π N ( x i ; μ 1 , σ 1 2 ) + ( 1 − π ) N ( x i ; μ 2 , σ 2 2 ) ) \log L(\theta)=\sum_{i}\log (\pi\mathcal N(x_i;\mu_1,\sigma^2_1)+(1-\pi)\mathcal N(x_i;\mu_2,\sigma^2_2)) logL(θ)=i∑log(πN(xi;μ1,σ12)+(1−π)N(xi;μ2,σ22))
其中参数 θ = ( π , μ 1 , σ 1 2 , μ 2 , σ 2 2 ) \theta=(\pi,\mu_1,\sigma^2_1,\mu_2,\sigma^2_2) θ=(π,μ1,σ12,μ2,σ22)。上式中存在"和的对数",若直接求导将会变得很麻烦。下节我们将会介绍EM算法来解决此类问题。
基本思想
概率模型有时既含有观测变量 (observable variable),又含有隐变量 (latent variable)。EM(Expectation-Maximization,期望最大算法)是一种迭代算法,用于含有隐变量的概率模型的极大似然估计或极大后验估计,是数据挖掘的十大经典算法之一。
假设现有一批独立同分布的样本
X = { x 1 , x 2 , ⋯ , x N } X=\{x_1,x_2,\cdots,x_N\} X={x1,x2,⋯,xN}
它们是由某个含有隐变量的概率分布 p ( x , z ∣ θ ) p(x,z|\theta) p(x,z∣θ) 生成。设样本对应的隐变量数据
Z = { z 1 , z 2 , ⋯ , z N } Z=\{z_1,z_2,\cdots,z_N\} Z={z1,z2,⋯,zN}
对于一个含有隐变量 Z Z Z 的概率模型,一般将 ( X , Z ) (X,Z) (X,Z) 称为完全数据 (complete-data),而观测数据 X X X 为不完全数据(incomplete-data)。
假设观测数据 X X X 概率密度函数是 p ( X ∣ θ ) p(X|\theta) p(X∣θ),其中 θ \theta θ是需要估计的模型参数,现尝试用极大似然估计法估计此概率分布的参数。为了便于讨论,此处假设 z z z 为连续型随机变量,则对数似然函数为
log L ( θ ) = ∑ i = 1 N log p ( x i ∣ θ ) = ∑ i = 1 N log ∫ z i p ( x i , z i ∣ θ ) d z i \log L(\theta)=\sum_{i=1}^N\log p(x_i|\theta)=\sum_{i=1}^N\log\int_{z_i}p(x_i,z_i|\theta)\mathrm dz_i logL(θ)=i=1∑Nlogp(xi∣θ)=i=1∑Nlog∫zip(xi,zi∣θ)dzi
Suppose you have a probability model with parameters θ \theta θ.
p ( x ∣ θ ) p(x|\theta) p(x∣θ) has two names. It can be called the probability of x x x (given θ \theta θ), or the likelihood of θ \theta θ (given that x x x was observed).
我们的目标是极大化观测数据 X X X 关于参数 θ \theta θ 的对数似然函数
θ ^ = arg max θ log L ( θ ) \hat\theta=\arg\max_{\theta}\log L(\theta) θ^=argθmaxlogL(θ)
显然,此时 log L ( θ ) \log L(\theta) logL(θ) 里含有未知的隐变量 z z z 以及求和项的对数,相比于不含隐变量的对数似然函数,该似然函数的极大值点较难求解,而 EM 算法则给出了一种迭代的方法来完成对 log L ( θ ) \log L(\theta) logL(θ) 的极大化。
注意:确定好含隐变量的模型后,即确定了联合概率密度函数 p ( x , z ∣ θ ) p(x,z|\theta) p(x,z∣θ) ,其中 θ \theta θ是需要估计的模型参数。为便于讨论,在此有必要说明下其他已知的概率函数。
联合概率密度函数
p ( x , z ∣ θ ) = f ( x , z ; θ ) p(x,z|\theta)=f(x,z;\theta) p(x,z∣θ)=f(x,z;θ)
观测变量 x x x 的概率密度函数
p ( x ∣ θ ) = ∫ z f ( x , z ; θ ) d z p(x|\theta)=\int_z f(x,z;\theta)\mathrm dz p(x∣θ)=∫zf(x,z;θ)dz
隐变量 z z z 的概率密度函数
p ( z ∣ θ ) = ∫ x f ( x , z ; θ ) d x p(z|\theta)=\int_x f(x,z;\theta)\mathrm dx p(z∣θ)=∫xf(x,z;θ)dx
条件概率密度函数
p ( x ∣ z , θ ) = p ( x , z ∣ θ ) p ( z ∣ θ ) = f ( x , z ; θ ) ∫ x f ( x , z ; θ ) d x p(x|z,\theta)=\frac{p(x,z|\theta)}{p(z|\theta)}=\frac{f(x,z;\theta)}{\int_x f(x,z;\theta)\mathrm dx} p(x∣z,θ)=p(z∣θ)p(x,z∣θ)=∫xf(x,z;θ)dxf(x,z;θ)
和
p ( z ∣ x , θ ) = p ( x , z ∣ θ ) p ( x ∣ θ ) = f ( x , z ; θ ) ∫ z f ( x , z ; θ ) d z p(z|x,\theta)=\frac{p(x,z|\theta)}{p(x|\theta)}=\frac{f(x,z;\theta)}{\int_z f(x,z;\theta)\mathrm dz} p(z∣x,θ)=p(x∣θ)p(x,z∣θ)=∫zf(x,z;θ)dzf(x,z;θ)
下面给出两种推导方法:一种借助 Jensen 不等式;一种使用 KL 散度。
首先使用 Jensen 不等式推导:使用含有隐变量的全概率公式
log p ( x i ∣ θ ) = log ∫ z i p ( x i , z i ∣ θ ) d z i = log ∫ z i q i ( z i ) p ( x i , z i ∣ θ ) q i ( z i ) d z i = log E z ( p ( x i , z i ∣ θ ) q i ( z i ) ) ⩾ E z ( log p ( x i , z i ∣ θ ) q i ( z i ) ) = ∫ z i q i ( z i ) log p ( x i , z i ∣ θ ) q i ( z i ) d z i \begin{aligned} \log p(x_i|\theta)&=\log\int_{z_i} p(x_i,z_i|\theta)\mathrm dz_i \\ &=\log\int_{z_i}q_i(z_i)\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\mathrm dz_i \\ &=\log\mathbb E_z\left(\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\right) \\ &\geqslant \mathbb E_z\left(\log\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\right) \\ &= \int_{z_i}q_i(z_i) \log\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\mathrm dz_i \end{aligned} logp(xi∣θ)=log∫zip(xi,zi∣θ)dzi=log∫ziqi(zi)qi(zi)p(xi,zi∣θ)dzi=logEz(qi(zi)p(xi,zi∣θ))⩾Ez(logqi(zi)p(xi,zi∣θ))=∫ziqi(zi)logqi(zi)p(xi,zi∣θ)dzi
其中 q i ( z i ) q_i(z_i) qi(zi) 是引入的第 i i i个样本隐变量 z i z_i zi 的任意概率密度函数(未知函数),其实 q q q 不管是任意函数,上式都成立。从后续推导得知,当 q i ( z i ) q_i(z_i) qi(zi) 是 z i z_i zi 的概率密度时,方便计算。
所以
log L ( θ ) = ∑ i = 1 N log p ( x i ∣ θ ) ⩾ B ( q , θ ) = ∑ i = 1 N ∫ z i q i ( z i ) log p ( x i , z i ∣ θ ) q i ( z i ) d z i \log L(\theta)=\sum_{i=1}^N\log p(x_i|\theta)\geqslant B(q,\theta)=\sum_{i=1}^N\int_{z_i}q_i(z_i) \log\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\mathrm dz_i logL(θ)=i=1∑Nlogp(xi∣θ)⩾B(q,θ)=i=1∑N∫ziqi(zi)logqi(zi)p(xi,zi∣θ)dzi
其中函数 B B B 为对数似然的下界函数。下界比较好求,所以我们要优化这个下界来使得似然函数最大。
假设第 t t t 次迭代时 θ \theta θ 的估计值是 θ ( t ) \theta^{(t)} θ(t),我们希望第 t + 1 t+1 t+1 次迭代时的 θ \theta θ 能使 log L ( θ ) \log L(\theta) logL(θ) 增大,即
log L ( θ ( t ) ) ⩽ log L ( θ ( t + 1 ) ) \log L(\theta^{(t)}) \leqslant \log L(\theta^{(t+1)}) logL(θ(t))⩽logL(θ(t+1))
可以分为两步实现:
-
首先,固定 θ = θ ( t ) \theta=\theta^{(t)} θ=θ(t) ,通过调整 q q q 函数使得 B ( q ( t ) , θ ) B(q^{(t)},\theta) B(q(t),θ) 在 θ ( t ) \theta^{(t)} θ(t) 处和 log L ( θ ( t ) ) \log L(\theta^{(t)}) logL(θ(t)) 相等;
B ( q ( t ) , θ ( t ) ) = log L ( θ ( t ) ) B(q^{(t)},\theta^{(t)})=\log L(\theta^{(t)}) B(q(t),θ(t))=logL(θ(t)) -
然后,固定 q q q,优化 θ ( t + 1 ) \theta^{(t+1)} θ(t+1) 取到下界函数 B ( q ( t ) , θ ) B(q^{(t)},\theta) B(q(t),θ) 的最大值。
θ ( t + 1 ) = arg max θ B ( q ( t ) , θ ) \theta^{(t+1)}=\arg\max_{\theta} B(q^{(t)},\theta) θ(t+1)=argθmaxB(q(t),θ)
所以
log L ( θ ( t + 1 ) ) ⩾ B ( q ( t ) , θ ( t + 1 ) ) ⩾ B ( q ( t ) , θ ( t ) ) = log L ( θ ( t ) ) \log L(\theta^{(t+1)})\geqslant B(q^{(t)},\theta^{(t+1)})\geqslant B(q^{(t)},\theta^{(t)})=\log L(\theta^{(t)}) logL(θ(t+1))⩾B(q(t),θ(t+1))⩾B(q(t),θ(t))=logL(θ(t))
因此,EM算法也可以看作一种坐标提升算法,首先固定一个值,对另外一个值求极值,不断重复直到收敛。
接下来,我们开始求解 q ( t ) q^{(t)} q(t) 。Jensen不等式中等号成立的条件是自变量是常数,即
p ( x i , z i ∣ θ ) q i ( z i ) = c \frac{p(x_i,z_i|\theta)}{q_i(z_i)}=c qi(zi)p(xi,zi∣θ)=c
由于假设 q i ( z i ) q_i(z_i) qi(zi)是 z i z_i zi 的概率密度函数,所以
p ( x i ∣ θ ) = ∫ z i p ( x i , z i ∣ θ ) d z i = ∫ z i c q i ( z i ) d z i = c p(x_i|\theta)=\int_{z_i}p(x_i,z_i|\theta)\mathrm dz_i=\int_{z_i} cq_i(z_i)\mathrm dz_i=c p(xi∣θ)=∫zip(xi,zi∣θ)dzi=∫zicqi(zi)dzi=c
于是
q i ( z i ) = p ( x i , z i ∣ θ ) c = p ( x i , z i ∣ θ ) p ( x i ∣ θ ) = p ( z i ∣ x i , θ ) q_i(z_i)=\frac{p(x_i,z_i|\theta)}{c}=\frac{p(x_i,z_i|\theta)}{p(x_i|\theta)}=p(z_i|x_i,\theta) qi(zi)=cp(xi,zi∣θ)=p(xi∣θ)p(xi,zi∣θ)=p(zi∣xi,θ)
可以看到,函数 q i ( z i ) q_i(z_i) qi(zi) 代表第 i i i 个数据是 z i z_i zi 的概率密度,是可以直接计算的。
最终,我们只要初始化或使用上一步已经固定的 θ ( t ) \theta^{(t)} θ(t),然后计算
θ ( t + 1 ) = arg max θ ∑ i = 1 N ∫ z i p ( z i ∣ x i , θ ( t ) ) log p ( x i , z i ∣ θ ) p ( z i ∣ x i , θ ( t ) ) d z i = arg max θ ∑ i = 1 N ∫ z i p ( z i ∣ x i , θ ( t ) ) log p ( x i , z i ∣ θ ) d z i = arg max θ ∑ i = 1 N E z i ∣ x i , θ ( t ) [ log p ( x i , z i ∣ θ ) ] = arg max θ Q ( θ , θ ( t ) ) \begin{aligned} \theta^{(t+1)}& = \arg\max_{\theta}\sum_{i=1}^N\int_{z_i}p(z_i|x_i,\theta^{(t)}) \log\frac{p(x_i,z_i|\theta)}{p(z_i|x_i,\theta^{(t)})}\mathrm dz_i \\ & = \arg\max_{\theta}\sum_{i=1}^N\int_{z_i}p(z_i|x_i,\theta^{(t)}) \log p(x_i,z_i|\theta)\mathrm dz_i \\ & = \arg\max_{\theta}\sum_{i=1}^N \mathbb E_{z_i|x_i,\theta^{(t)}}[\log p(x_i,z_i|\theta)] \\ & = \arg\max_{\theta} Q(\theta,\theta^{(t)}) \end{aligned} θ(t+1)=argθmaxi=1∑N∫zip(zi∣xi,θ(t))logp(zi∣xi,θ(t))p(xi,zi∣θ)dzi=argθmaxi=1∑N∫zip(zi∣xi,θ(t))logp(xi,zi∣θ)dzi=argθmaxi=1∑NEzi∣xi,θ(t)[logp(xi,zi∣θ)]=argθmaxQ(θ,θ(t))
接下来使用 KL 散度推导:使用含有隐变量的条件概率
log p ( x i ∣ θ ) = log p ( x i , z i ∣ θ ) p ( z i ∣ x i , θ ) = ∫ z i q i ( z i ) log p ( x i , z i ∣ θ ) p ( z i ∣ x i , θ ) ⋅ q i ( z i ) q i ( z i ) d z i = ∫ z i q i ( z i ) log p ( x i , z i ∣ θ ) q i ( z i ) d z i + ∫ z i q i ( z i ) log q i ( z i ) p ( z i ∣ x i , θ ) d z i = B ( q i , θ ) + K L ( q i ∥ p i ) \begin{aligned} \log p(x_i|\theta)&=\log\frac{p(x_i,z_i|\theta)}{p(z_i|x_i,\theta)} \\ &=\int_{z_i}q_i(z_i)\log\frac{p(x_i,z_i|\theta)}{p(z_i|x_i,\theta)}\cdot\frac{q_i(z_i)}{q_i(z_i)}\mathrm dz_i \\ &= \int_{z_i}q_i(z_i) \log\frac{p(x_i,z_i|\theta)}{q_i(z_i)}\mathrm dz_i + \int_{z_i}q_i(z_i) \log\frac{q_i(z_i)}{p(z_i|x_i,\theta)}\mathrm dz_i \\ &=B(q_i,\theta)+KL(q_i\|p_i) \end{aligned} logp(xi∣θ)=logp(zi∣xi,θ)p(xi,zi∣θ)=∫ziqi(zi)logp(zi∣xi,θ)p(xi,zi∣θ)⋅qi(zi)qi(zi)dzi=∫ziqi(zi)logqi(zi)p(xi,zi∣θ)dzi+∫ziqi(zi)logp(zi∣xi,θ)qi(zi)dzi=B(qi,θ)+KL(qi∥pi)
同样 q i ( z i ) q_i(z_i) qi(zi) 是引入的关于 z i z_i zi 的任意概率密度函数(未知函数),函数 B ( q i , θ ) B(q_i,\theta) B(qi,θ) 表示对数似然的一个下界,散度 K L ( q i ∥ p i ) KL(q_i\|p_i) KL(qi∥pi)描述了下界与对数似然的差距。
同样为了保证
log L ( θ ( t ) ) ⩽ log L ( θ ( t + 1 ) ) \log L(\theta^{(t)}) \leqslant \log L(\theta^{(t+1)}) logL(θ(t))⩽logL(θ(t+1))
分为两步实现:
-
首先,固定 θ = θ ( t ) \theta=\theta^{(t)} θ=θ(t) ,通过调整 q q q 函数使得 B ( q ( t ) , θ ) B(q^{(t)},\theta) B(q(t),θ) 在 θ ( t ) \theta^{(t)} θ(t) 处和 log L ( θ ( t ) ) \log L(\theta^{(t)}) logL(θ(t)) 相等,即 K L ( q i ∥ p i ) = 0 KL(q_i\|p_i)=0 KL(qi∥pi)=0,于是
q i ( z i ) = p ( z i ∣ x i , θ ( t ) ) q_i(z_i)=p(z_i|x_i,\theta^{(t)}) qi(zi)=p(zi∣xi,θ(t)) -
然后,固定 q q q,优化 θ ( t + 1 ) \theta^{(t+1)} θ(t+1) 取到下界函数 B ( q ( t ) , θ ) B(q^{(t)},\theta) B(q(t),θ) 的最大值。
θ ( t + 1 ) = arg max θ B ( q ( t ) , θ ) \theta^{(t+1)}=\arg\max_{\theta} B(q^{(t)},\theta) θ(t+1)=argθmaxB(q(t),θ)
算法流程
输入:观测数据 X X X,联合分布 p ( x , z ; θ ) p(x,z;\theta) p(x,z;θ),条件分布 P ( z ∣ x , θ ) P(z|x,\theta) P(z∣x,θ)
输出:模型参数 θ \theta θ
EM算法通过引入隐含变量,使用极大似然估计(MLE)进行迭代求解参数。每次迭代由两步组成:
-
E-step:求期望 (expectation)。以参数的初始值或上一次迭代的模型参数 θ ( t ) \theta^{(t)} θ(t) 来计算隐变量后验概率 p ( z i ∣ x i , θ ( t ) ) p(z_i|x_i,\theta^{(t)}) p(zi∣xi,θ(t)) ,并计算期望(expectation)
Q ( θ , θ ( t ) ) = ∑ i = 1 N ∫ z i p ( z i ∣ x i , θ ( t ) ) log p ( x i , z i ∣ θ ) d z i Q(\theta,\theta^{(t)})=\sum_{i=1}^N\int_{z_i}p(z_i|x_i,\theta^{(t)}) \log p(x_i,z_i|\theta)\mathrm dz_i Q(θ,θ(t))=i=1∑N∫zip(zi∣xi,θ(t))logp(xi,zi∣θ)dzi -
M-step: 求极大 (maximization),极大化E步中的期望值,来确定 t + 1 t+1 t+1 次迭代的参数估计值
θ ( t + 1 ) = arg max θ Q ( θ , θ ( t ) ) \theta^{(t+1)}=\arg\max_{\theta} Q(\theta,\theta^{(t)}) θ(t+1)=argθmaxQ(θ,θ(t))
依次迭代,直至收敛到局部最优解。
高斯混合模型
基础模型
高斯混合模型 (Gaussian Mixture Model, GMM) 数据可以看作是从 K K K个高斯分布中生成出来的,每个高斯分布称为一个组件 (Component)。
引入隐变量 z ∈ { 1 , 2 , ⋯ , K } z\in\{1,2,\cdots,K\} z∈{1,2,⋯,K},表示对应的样本 x x x 属于哪一个高斯分布,这个变量是一个离散的随机变量:
P ( z = k ) = π k s.t. ∑ k = 1 K π k = 1 \mathbb P(z=k)=\pi_k \\ \text{s.t. } \sum_{k=1}^K\pi_k=1 P(z=k)=πks.t. k=1∑Kπk=1
可将 π k \pi_k πk 视为选择第 k k k 高斯分布的先验概率,而对应的第 k k k 个高斯分布的样本概率
p ( x ∣ z = k ) = N ( x ; μ k , Σ k ) p(x|z=k)=\mathcal N(x;\mu_k,\Sigma_k) p(x∣z=k)=N(x;μk,Σk)
于是高斯混合模型
p M ( x ) = ∑ k = 1 K π k N ( x ; μ k , Σ k ) p_M(x)=\sum_{k=1}^K\pi_k\mathcal N(x;\mu_k,\Sigma_k) pM(x)=k=1∑KπkN(x;μk,Σk)
其中 0 ⩽ π k ⩽ 1 0\leqslant \pi_k\leqslant 1 0⩽πk⩽1为混合系数(mixing coefficients)。
高斯混合模型的参数估计是EM算法的一个重要应用,隐马尔科夫模型的非监督学习也是EM算法的一个重要应用。
EM算法
高斯混合模型的极大似然估计
θ ^ = arg max θ ∑ i = 1 N log ∑ k = 1 K π k N ( x i ; μ k , Σ k ) \hat\theta=\arg\max_{\theta} \sum_{i=1}^N\log\sum_{k=1}^K\pi_k \mathcal N(x_i;\mu_k,\Sigma_k) θ^=argθmaxi=1∑Nlogk=1∑KπkN(xi;μk,Σk)
其中参数 θ k = ( π k , μ k , Σ k ) \theta_k=(\pi_k,\mu_k,\Sigma_k) θk=(πk,μk,Σk),使用EM算法估计GMM的参数 θ \theta θ。
依照当前模型参数,计算隐变量后验概率:由贝叶斯公式知道
P ( z i = k ∣ x i ) = P ( z i = k ) p ( x i ∣ z i = k ) p ( x i ) = π k N ( x i ; μ k , Σ k ) ∑ k = 1 K π k N ( x i ; μ k , Σ k ) = γ i k \begin{aligned} P(z_i=k|x_i)&=\frac{P(z_i=k)p(x_i|z_i=k)}{p(x_i)} \\ &=\frac{\pi_k\mathcal N(x_i;\mu_k,\Sigma_k)}{\sum_{k=1}^K\pi_k\mathcal N(x_i;\mu_k,\Sigma_k) } \\ &=\gamma_{ik} \end{aligned} P(zi=k∣xi)=p(xi)P(zi=k)p(xi∣zi=k)=∑k=1KπkN(xi;μk,Σk)πkN(xi;μk,Σk)=γik
令 γ i k \gamma_{ik} γik表示第 i i i个样本属于第 k k k个高斯分布的概率。
E-step:确定Q函数
Q ( θ , θ ( t ) ) = ∑ i = 1 N ∑ k = 1 K p ( z i = k ∣ x i , μ ( t ) , Σ ( t ) ) log p ( x i , z i = k ∣ μ , Σ ) = ∑ i = 1 N ∑ k = 1 K γ i k log π k N ( x ; μ k , Σ k ) = ∑ i = 1 N ∑ k = 1 K γ i k ( log π k + log N ( x ; μ k , Σ k ) ) \begin{aligned} Q(\theta,\theta^{(t)})&=\sum_{i=1}^N\sum_{k=1}^Kp(z_i=k|x_i,\mu^{(t)},\Sigma^{(t)}) \log p(x_i,z_i=k|\mu,\Sigma) \\ &=\sum_{i=1}^N\sum_{k=1}^K\gamma_{ik}\log\pi_k\mathcal N(x;\mu_k,\Sigma_k) \\ &=\sum_{i=1}^N\sum_{k=1}^K\gamma_{ik}(\log\pi_k+ \log\mathcal N(x;\mu_k,\Sigma_k) ) \end{aligned} Q(θ,θ(t))=i=1∑Nk=1∑Kp(zi=k∣xi,μ(t),Σ(t))logp(xi,zi=k∣μ,Σ)=i=1∑Nk=1∑KγiklogπkN(x;μk,Σk)=i=1∑Nk=1∑Kγik(logπk+logN(x;μk,Σk))
M-step:求Q函数的极大值
上面已获得的 Q ( θ , θ ( t ) ) Q(\theta,\theta^{(t)}) Q(θ,θ(t))分别对 μ k , Σ k \mu_k,\Sigma_k μk,Σk求导并设为0。得到
μ k ( t + 1 ) = ∑ i = 1 N γ i k x i ∑ i = 1 N γ i k Σ k ( t + 1 ) = ∑ i = 1 N γ i k ( x i − μ k ( t + 1 ) ) ( x i − μ k ( t + 1 ) ) T ∑ i = 1 N γ i k \mu_k^{(t+1)}=\frac{\sum_{i=1}^N\gamma_{ik}x_i}{\sum_{i=1}^N\gamma_{ik}} \\ \Sigma_k^{(t+1)}=\frac{\sum_{i=1}^N\gamma_{ik}(x_i-\mu_k^{(t+1)}) (x_i-\mu_k^{(t+1)})^T }{\sum_{i=1}^N\gamma_{ik}} μk(t+1)=∑i=1Nγik∑i=1NγikxiΣk(t+1)=∑i=1Nγik∑i=1Nγik(xi−μk(t+1))(xi−μk(t+1))T
可以看到第 k k k个高斯分布的 μ k , Σ k \mu_k,\Sigma_k μk,Σk 是所有样本的加权平均,其中每个样本的权重为该样本属于第 k k k个高斯分布的后验概率 γ i k \gamma_{ik} γik。
对于混合系数 π k \pi_k πk,因为有限制条件,使用拉格朗日乘子法可求得
π k ( t + 1 ) = 1 N ∑ i = 1 N γ i k \pi_k^{(t+1)}=\frac{1}{N}\sum_{i=1}^N\gamma_{ik} πk(t+1)=N1i=1∑Nγik
即第 k k k个高斯分布的混合系数是属于 k k k的样本的平均后验概率,由此运用EM算法能大大简化高斯混合模型的参数估计过程,在中间步只需计算 γ i k \gamma_{ik} γik就行了。
高斯混合模型的算法流程如下图所示:
高斯混合聚类
高斯混合聚类假设每个类簇中的样本都服从一个多维高斯分布,那么空间中的样本可以看作由 K K K个多维高斯分布混合而成。
引入隐变量 z z z 标记簇类别,这样就可以使用高斯混合模型
p M ( x ) = ∑ k = 1 K π k N ( x ; μ k , Σ k ) p_M(x)=\sum_{k=1}^K\pi_k\mathcal N(x;\mu_k,\Sigma_k) pM(x)=k=1∑KπkN(x;μk,Σk)
使用EM算法迭代求解。
相比于K-means更具一般性,能形成各种不同大小和形状的簇。K-means可视为高斯混合聚类中每个样本仅指派给一个混合成分的特例
γ i k = { 1 , if k = arg min k ∥ x i − μ k ∥ 2 0 , otherwise \gamma_{ik}=\begin{cases} 1, & \text{if } k=\arg\min_k\|x_i-\mu_k\|^2\\ 0, & \text{otherwise} \end{cases} γik={1,0,if k=argmink∥xi−μk∥2otherwise
且各混合成分协方差相等,均为对角矩阵 σ 2 I \sigma^2 I σ2I。
附录
Jensen 不等式
若 f f f 是凸函数(convex function),对任意的 λ ∈ [ 0 , 1 ] \lambda\in [0,1] λ∈[0,1],下式恒成立
f ( λ x 1 + ( 1 − λ ) x 2 ) ⩽ λ f ( x 1 ) + ( 1 − λ ) f ( x 2 ) f(\lambda x_1+(1-\lambda)x_2)\leqslant \lambda f(x_1)+(1-\lambda)f(x_2) f(λx1+(1−λ)x2)⩽λf(x1)+(1−λ)f(x2)
Jensen’s inequality就是上式的推广,设 f ( x ) f(x) f(x) 为凸函数, λ i ∈ [ 0 , 1 ] , ∑ i λ i = 1 \lambda_i\in[0,1],\ \sum_i\lambda_i=1 λi∈[0,1], ∑iλi=1 ,则
f ( ∑ i λ i x i ) ⩽ ∑ i λ i f ( x i ) f(\sum_i\lambda_ix_i)\leqslant \sum_i\lambda_if(x_i) f(i∑λixi)⩽i∑λif(xi)
若将 λ i \lambda_i λi 视为一个概率分布,则可表示为期望值的形式
f ( E [ x ] ) ⩽ E [ f ( x ) ] f(\mathbb E[x])\leqslant\mathbb E[f(x)] f(E[x])⩽E[f(x)]
显然,如果 f f f 是凹函数(concave function),则将不等号反向。