机器学习——EM算法


以下内容均为个人理解,如有错误,欢迎指出


什么是EM算法

EM算法是一类通过迭代进行极大似然估计的优化算法,用于对包含隐变量或缺失数据的概率模型进行参数估计。


什么是隐变量

可以表示缺失数据,或概率模型中任何无法直接观测的随机变量,隐变量往往会影响观测变量的取值,例如现在有一群人的身高,但是我们并不知道每个人的性别,性别就是一个隐变量

什么是极大似然估计

  • 概念:我们有一个概率密度(质量)函数,其中含有未知数,极大似然估计就是利用已知的样本结果,反推最有可能(最大概率)导致这样结果的未知数的值。
  • 可以使用最大似然估计的条件:

1、n次实验互相独立,即样本之间互相独立
2、已知其概率密度函数(连续分布)或概率质量函数(离散分布)

  • 举个例子,对于某个事件,我们已知其概率质量函数(离散分布):P(Y1)、P(Y2)… P(Yn),即给定一个输入,输出为Y1、Y2… Yn的概率,但是这些概率质量函数中有一个未知数,记为W,那么如何确定W的值呢?我们给定n个输入,进行n次实验,n次实验之间相互独立,获得n个输出,记为Y,极大似然估计相信进行n次实验,最终结果为Y的概率最大,由于n次实验之间互相独立,记第i次实验结果对应的概率质量函数为Zi,那么结果为Y的概率为 ∏ i = 1 n \prod_{i=1}^n i=1nZi,此时就是一个最优化问题,即确定W,使 ∏ i = 1 n \prod_{i=1}^n i=1nZi有最大值,最优化那一套理论便可以派上用场,由于乘法求导比较麻烦,所以我们常常加一个log,即log( ∏ i = 1 n \prod_{i=1}^n i=1nZi),此时可变为加法求导。

为什么需要EM算法

以下为个人理解,如有错误,欢迎指出

给定一个南方人身高的观测样本 x 1 , x 2 , x 3 , . . . . x n x_1,x_2,x_3,....x_n x1,x2,x3,....xn,假定身高服从正态分布,我们不知道正态分布的期望与方差,这两个参数记为 θ \theta θ,现在我们想依据这批身高数据求得正态分布的两个参数,从而求得南方人身高的概率密度函数。

我们知道性别可以间接影响身高,例如一个人是女性,那么大概率其身高在1米6左右,如果是男性,其身高大概率在1米7左右,如果我们直接使用极大似然估计最大化对数似然 m a x ( ∑ i = 1 n l o g p ( x i ∣ θ ) ) max(\sum_{i=1}^nlogp(x_i|\theta)) max(i=1nlogp(xiθ)),假设这批样本都是女性,那么求得的概率密度函数将显示大部分南方人的身高都是1米6(男孩子哭死),显然,在有隐变量的情况下,直接使用极大似然估计求出的参数 θ \theta θ将具有很大的偏差,而EM算法将在考虑隐变量的情况下最大化对数似然函数,则有
m a x ( ∑ i = 1 n log ⁡ p ( x i ∣ θ ) ) = m a x ( ∑ i = 1 n log ⁡ ∑ j = 1 m p ( x i , z j ∣ θ ) ) ( 隐 变 量 是 离 散 随 机 变 量 ) m a x ( ∑ i = 1 n log ⁡ p ( x i ∣ θ ) ) = m a x ( ∑ i = 1 n log ⁡ ∫ a b p ( x i , z ∣ θ ) d z ) ( 隐 变 量 是 连 续 随 机 变 量 ) \begin{aligned} max(\sum_{i=1}^n\log p(x_i|\theta))=&max(\sum_{i=1}^n\log \sum_{j=1}^mp(x_i,z_j|\theta))(隐变量是离散随机变量)\\ max(\sum_{i=1}^n\log p(x_i|\theta))=&max(\sum_{i=1}^n\log \int_a^bp(x_i,z|\theta)dz)(隐变量是连续随机变量) \end{aligned} max(i=1nlogp(xiθ))=max(i=1nlogp(xiθ))=max(i=1nlogj=1mp(xi,zjθ))max(i=1nlogabp(xi,zθ)dz)
这看起来非常神奇,毕竟我们都不知道样本集中的样本的性别,但是EM算法是有理论保障的,本文将假设隐变量为离散随机变量


EM算法的推导

预备知识

Jensen不等式
x是一个随机变量,f(x)是一个概率密度函数,如果它是一个凸函数,则有
E [ f ( x ) ] ≥ f [ E ( x ) ] E[f(x)]\geq f[E(x)] E[f(x)]f[E(x)]

若f(x)是一个凹函数,则有
E [ f ( x ) ] ≤ f [ E ( x ) ] E[f(x)]\leq f[E(x)] E[f(x)]f[E(x)]

当且仅当x是一个常数时,即随机变量取值固定,等号成立


推导过程

符号表

符号名含义
Q ( z i j ) Q(z_{ij}) Q(zij)表示样本i属于隐藏变量 z j z_j zj的概率
m样本数
n隐变量数
θ \theta θ未知参数
p ( x i , z j ; θ ) p(x_i,z_j;\theta) p(xi,zj;θ)表示隐变量取值为 z j z_j zj且样本 x i x_i xi出现的概率(带未知参数的表达式)
y j y_j yj y j = p ( x i , z j ; θ ) Q ( z i j ) y_j=\frac{p(x_i,z_j;\theta)}{Q(z_{ij})} yj=Q(zij)p(xi,zj;θ),为离散随机变量

我们的想法是构造出Jenson不等式,对于对数似然 ∑ i = 1 m log ⁡ ∑ j = 1 n p ( x i , z j ; θ ) \sum_{i=1}^m\log \sum_{j=1}^np(x_i,z_j;\theta) i=1mlogj=1np(xi,zj;θ),我们取其中的一项 log ⁡ ∑ j = 1 n p ( x i , z j ; θ ) \log \sum_{j=1}^np(x_i,z_j;\theta) logj=1np(xi,zj;θ)进行讨论,对其进行变形,则有
log ⁡ ∑ j = 1 n p ( x i , z j ; θ ) = log ⁡ ∑ j = 1 n ( Q ( z i j ) p ( x i , z j ; θ ) Q ( z i j ) ) ( 式 1.1 ) \log\sum_{j=1}^np(x_i,z_j;\theta)=\log \sum_{j=1}^n(Q(z_{ij})\frac{p(x_i,z_j;\theta)}{Q(z_{ij})})(式1.1) logj=1np(xi,zj;θ)=logj=1n(Q(zij)Q(zij)p(xi,zj;θ))1.1
若设 p ( Y = y j ) = Q ( z i j ) p(Y=y_j)=Q(z_{ij}) p(Y=yj)=Q(zij),式1.1可看成 ∑ j = 1 n p ( Y = y j ) y j \sum_{j=1}^np(Y=y_j)y_j j=1np(Y=yj)yj,这是离散随机变量 y j y_j yj的期望,假设有一个对数函数,则式1.1可看成 f ( E ( Y ) ) f(E(Y)) f(E(Y)),而对数函数是凹函数,所以 E ( f ( Y ) ) = ∑ i = j n [ l o g ( p ( x i , z j ; θ ) Q ( z i j ) ) ] Q ( z i j ) E(f(Y))=\sum_{i=j}^n[log(\frac{p(x_i,z_j;\theta)}{Q(z_{ij})})]Q(z_{ij}) E(f(Y))=i=jn[log(Q(zij)p(xi,zj;θ))]Q(zij),由于对数函数是凹函数,则有 f ( E ( Y ) ) = log ⁡ ∑ j = 1 n ( Q ( z i j ) p ( x i , z j ; θ ) Q ( z i j ) ) ≥ E ( f ( Y ) ) = ∑ j = 1 n [ l o g ( p ( x i , z j ; θ ) Q ( z i j ) ) ] Q ( z i j ) f(E(Y))=\log \sum_{j=1}^n(Q(z_{ij})\frac{p(x_i,z_j;\theta)}{Q(z_{ij})})\geq E(f(Y))=\sum_{j=1}^n[log(\frac{p(x_i,z_j;\theta)}{Q(z_{ij})})]Q(z_{ij}) f(E(Y))=logj=1n(Q(zij)Q(zij)p(xi,zj;θ))E(f(Y))=j=1n[log(Q(zij)p(xi,zj;θ))]Q(zij)
等号成立的条件是 y j   ( 1 ≤ j ≤ n ) y_j\ (1 \leq j \leq n) yj (1jn)均为常数,设该常数为c,则有
y 1 = p ( x i , z 1 ; θ ) Q ( z i 1 ) = c y 2 = p ( x i , z 2 ; θ ) Q ( z i 2 ) = c . . . . . y n = p ( x i , z 1 ; θ ) Q ( z i n ) = c \begin{aligned} y_1&=\frac{p(x_i,z_1;\theta)}{Q(z_{i1})}=c\\ y_2&=\frac{p(x_i,z_2;\theta)}{Q(z_{i2})}=c\\ &.....\\ y_n&=\frac{p(x_i,z_1;\theta)}{Q(z_{in})}=c \end{aligned} y1y2yn=Q(zi1)p(xi,z1;θ)=c=Q(zi2)p(xi,z2;θ)=c.....=Q(zin)p(xi,z1;θ)=c
则有
p ( x i , z 1 ; θ ) c = Q ( z i 1 ) p ( x i , z 2 ; θ ) c = Q ( z i 2 ) . . . . . p ( x i , z n ; θ ) c = Q ( z i n ) \begin{aligned} \frac{p(x_i,z_1;\theta)}{c}&=Q(z_{i1})\\ \frac{p(x_i,z_2;\theta)}{c}&=Q(z_{i2})\\ .....\\ \frac{p(x_i,z_n;\theta)}{c}&=Q(z_{in}) \end{aligned} cp(xi,z1;θ)cp(xi,z2;θ).....cp(xi,zn;θ)=Q(zi1)=Q(zi2)=Q(zin)
由于一个样本必然属于某个隐变量,则有 ∑ i = 1 n Q ( z i j ) = 1 \sum_{i=1}^nQ(z_{ij})=1 i=1nQ(zij)=1,将上述式子相加即为
c = ∑ j = 1 n p ( x i , z j ; θ ) = ∑ j = 1 n p ( x i ; θ ) p ( z j ∣ x i ; θ ) = p ( x i ; θ ) ∑ j = 1 n p ( z j ∣ x i ; θ ) = p ( x i ; θ ) c=\sum_{j=1}^np(x_i,z_j;\theta)=\sum_{j=1}^np(x_i;\theta)p(z_j|x_i;\theta)=p(x_i;\theta)\sum_{j=1}^np(z_j|x_i;\theta)=p(x_i;\theta) c=j=1np(xi,zj;θ)=j=1np(xi;θ)p(zjxi;θ)=p(xi;θ)j=1np(zjxi;θ)=p(xi;θ)
所以有
Q ( z i j ) = p ( x i , z j ; θ ) p ( x i ; θ ) = p ( z j ∣ x i ; θ ) Q(z_{ij})=\frac{p(x_i,z_j;\theta)}{p(x_i;\theta)}=p(z_j|x_i;\theta) Q(zij)=p(xi;θ)p(xi,zj;θ)=p(zjxi;θ)
即已知样本i,该样本属于隐变量 z j z_j zj的概率,此时等号成立,即
log ⁡ ∑ j = 1 n p ( x i , z j ; θ ) = ∑ j = 1 n [ log ⁡ ( p ( x i , z j ; θ ) Q ( z i j ) ) ] Q ( z i j ) = ∑ j = 1 n [ log ⁡ p ( x i , z j ; θ ) p ( z j ∣ x i ; θ ) ] p ( z j ∣ x i ; θ ) \log \sum_{j=1}^n p(x_i,z_j;\theta)=\sum_{j=1}^n[\log(\frac{p(x_i,z_j;\theta)}{Q(z_{ij})})]Q(z_{ij})=\sum_{j=1}^n[\log \frac{p(x_i,z_j;\theta)}{p(z_j|x_i;\theta)}]p(z_j|x_i;\theta) logj=1np(xi,zj;θ)=j=1n[log(Q(zij)p(xi,zj;θ))]Q(zij)=j=1n[logp(zjxi;θ)p(xi,zj;θ)]p(zjxi;θ)
则有
∑ i = 1 m log ⁡ ∑ j = 1 n p ( x i , z j ; θ ) = ∑ i = 1 m ∑ j = 1 n [ log ⁡ p ( x i , z j ; θ ) p ( z j ∣ x i ; θ ) ] p ( z j ∣ x i ; θ ) ( 式 1.2 ) \sum_{i=1}^m\log \sum_{j=1}^n p(x_i,z_j;\theta)=\sum_{i=1}^m\sum_{j=1}^n[\log \frac{p(x_i,z_j;\theta)}{p(z_j|x_i;\theta)}]p(z_j|x_i;\theta)(式1.2) i=1mlogj=1np(xi,zj;θ)=i=1mj=1n[logp(zjxi;θ)p(xi,zj;θ)]p(zjxi;θ)1.2

因此最大化对数似然等价于最大化式1.2的右半部分(以下简称式1.2),为什么要使用Jenson不等式进行变化呢?这是由于对数似然在对数函数内部有一个和式,加上概率密度函数一般非常复杂,例如正态分布,求导的结果将会非常复杂,导致求解最大值异常困难
虽然式1.2看起来更加复杂,但是由于EM算法的存在,求其(局部)最大的过程将较为简洁


EM算法的过程

EM算法是一个迭代过程,假设我们的概率密度函数的一组未知参数为 θ 1 , θ 2 , . . . , θ N \theta_1,\theta_2,...,\theta_N θ1,θ2,...,θN,记为 θ \theta θ,我们首先随机初始化这组参数,假设有k个样本,隐藏变量为 z 1 , z 2 , . . . . , z n z_1,z_2,....,z_n z1,z2,....,zn,设第t次迭代,求得的参数取值为 θ 1 t , θ 2 t , . . . , θ N t \theta_1^t,\theta_2^t,...,\theta_N^t θ1t,θ2t,...,θNt,记为 θ t \theta^t θt,则第t+1次迭代的参数取值 θ 1 t + 1 , θ 2 t + 1 , . . . . . , θ N t + 1 \theta_1^{t+1},\theta_2^{t+1},.....,\theta_N^{t+1} θ1t+1,θ2t+1,.....,θNt+1(记为 θ t + 1 \theta^{t+1} θt+1)通过下列步骤求得

  • E步:根据参数 θ t , \theta^t, θt,计算每个样本属于 z j   ( 1 ≤ j ≤ n ) z_j \ (1\leq j \leq n) zj (1jn)的概率,即 p ( z j ∣ x i ; θ t ) p(z_j|x_i;\theta^t) p(zjxi;θt),计算下式:
    ∑ i = 1 k ∑ j = 1 n p ( z j ∣ x i ; θ t ) l o g p ( x i , z j ; θ ) p ( z j ∣ x i ; θ t ) ( 式 1.3 ) \sum_{i=1}^k\sum_{j=1}^n p(z_j|x_i;\theta^t) log\frac{p(x_i,z_j;\theta)}{p(z_j|x_i;\theta^t)} (式1.3) i=1kj=1np(zjxi;θt)logp(zjxi;θt)p(xi,zj;θ)1.3其中, p ( x i , z j ; θ ) p(x_i,z_j;\theta) p(xi,zj;θ)表示隐变量取值为 z j z_j zj且样本 x i x_i xi出现的概率,注意到参数是 θ \theta θ,是变量,式1.3可以看成是复合随机变量的期望。
  • M步:根据E步的结果,接下来我们将最大化式1.3,此时 θ \theta θ的取值即为 θ t + 1 \theta^{t+1} θt+1

重复E步与M步,直到参数不再变化为止,相比于直接对对数似然求导,式1.3求导的结果将非常简洁。

总结一下为什么需要使用Jenson不等式对含有隐变量的对数似然进行转换。
因为含有隐变量的对数似然求导非常困难,这导致求解最大化过程非常困难,因此,我们使用Jenson不等式对其进行转换,接着用EM算法进行优化,此时求解过程将变得更加简洁容易


EM算法为什么有效

我们只需要证明EM算法的每次迭代求得的 θ \theta θ都会使式1.2增大即可,即
∑ i = 1 n log ⁡ p ( x i ; θ t + 1 ) = ∑ i = 1 k ∑ j = 1 n [ log ⁡ p ( x i , z j ; θ t + 1 ) p ( z j ∣ x i ; θ t + 1 ) ] p ( z j ∣ x i ; θ t + 1 ) ≥ ∑ i = 1 k ∑ j = 1 n [ log ⁡ p ( x i , z j ; θ t ) p ( z j ∣ x i ; θ t ) ] p ( z j ∣ x i ; θ t ) = ∑ i = 1 n log ⁡ p ( x i ; θ t ) ( 式 1.5 ) \sum_{i=1}^n \log p(x_i;\theta^{t+1})=\sum_{i=1}^k\sum_{j=1}^n[\log \frac{p(x_i,z_j;\theta^{t+1})}{p(z_j|x_i;\theta^{t+1})}]p(z_j|x_i;\theta^{t+1})\geq\sum_{i=1}^k\sum_{j=1}^n[\log \frac{p(x_i,z_j;\theta^t)}{p(z_j|x_i;\theta^t)}]p(z_j|x_i;\theta^t)=\sum_{i=1}^n \log p(x_i;\theta^t)(式1.5) i=1nlogp(xi;θt+1)=i=1kj=1n[logp(zjxi;θt+1)p(xi,zj;θt+1)]p(zjxi;θt+1)i=1kj=1n[logp(zjxi;θt)p(xi,zj;θt)]p(zjxi;θt)=i=1nlogp(xi;θt)1.5
其中 θ t 、 θ t + 1 \theta^t、\theta^{t+1} θtθt+1分别表示第t、t+1轮变量 θ \theta θ的取值,等号成立时表示位于(局部)最大

依据条件概率,我们有
p ( x ; θ ) = p ( x , z ; θ ) p ( z ∣ x ; θ ) p(x;\theta)=\frac{p(x,z;\theta)}{p(z|x;\theta)} p(x;θ)=p(zx;θ)p(x,z;θ)
取对数可得
log ⁡ p ( x ; θ ) = log ⁡ p ( x , z ∣ θ ) − log ⁡ p ( z ∣ x ; θ ) \log p(x;\theta)=\log p(x,z|\theta)-\log p(z|x;\theta) logp(x;θ)=logp(x,zθ)logp(zx;θ)

Q ( θ , θ t ) = ∑ i = 1 n ∑ j = 1 m log ⁡ p ( x i , z i j ; θ ) p ( z i j ∣ x i ; θ t ) H ( θ , θ t ) = ∑ i = 1 n ∑ j = 1 m log ⁡ p ( z i j ∣ x i ; θ ) p ( z i j ∣ x i ; θ t ) \begin{aligned} Q(\theta,\theta^t)=&\sum_{i=1}^n\sum_{j=1}^m\log p(x_i,z_{ij};\theta)p(z_{ij}|x_i;\theta^t)\\ H(\theta,\theta^t)=&\sum_{i=1}^n\sum_{j=1}^m\log p(z_{ij}|x_i;\theta)p(z_{ij}|x_i;\theta^t) \end{aligned} Q(θ,θt)=H(θ,θt)=i=1nj=1mlogp(xi,zij;θ)p(zijxi;θt)i=1nj=1mlogp(zijxi;θ)p(zijxi;θt)
所以有
Q ( θ , θ t ) − H ( θ , θ t ) = ∑ i = 1 n ∑ j = 1 m [ log ⁡ p ( x i , z i j ; θ ) − log ⁡ p ( z i j ∣ x i ; θ ) ] p ( z i j ∣ x i ; θ t ) = ∑ i = 1 n ∑ j = 1 m log ⁡ p ( x i , z i j ; θ ) p ( z i j ∣ x i ; θ ) p ( z i j ∣ x i ; θ t ) = ∑ i = 1 n ∑ j = 1 m log ⁡ p ( x i ; θ ) p ( z i j ∣ x i ; θ t ) ( 条 件 概 率 的 公 式 ) = ∑ i = 1 n log ⁡ p ( x i ; θ ) ∑ j = 1 m p ( z i j ∣ x i ; θ t ) = ∑ i = 1 n log ⁡ p ( x i ; θ ) \begin{aligned} Q(\theta,\theta^t)-H(\theta,\theta^t)&=\sum_{i=1}^n\sum_{j=1}^m[\log p(x_i,z_{ij};\theta)-\log p(z_{ij}|x_i;\theta)]p(z_{ij}|x_i;\theta^t)\\ &=\sum_{i=1}^n\sum_{j=1}^m \log \frac{p(x_i,z_{ij};\theta)}{p(z_{ij}|x_i;\theta)} p(z_{ij}|x_i;\theta^t)\\ &=\sum_{i=1}^n\sum_{j=1}^m \log p(x_i;\theta) p(z_{ij}|x_i;\theta^t)(条件概率的公式)\\ &=\sum_{i=1}^n \log p(x_i;\theta) \sum_{j=1}^m p(z_{ij}|x_i;\theta^t)\\ &=\sum_{i=1}^n \log p(x_i;\theta) \end{aligned} Q(θ,θt)H(θ,θt)=i=1nj=1m[logp(xi,zij;θ)logp(zijxi;θ)]p(zijxi;θt)=i=1nj=1mlogp(zijxi;θ)p(xi,zij;θ)p(zijxi;θt)=i=1nj=1mlogp(xi;θ)p(zijxi;θt)()=i=1nlogp(xi;θ)j=1mp(zijxi;θt)=i=1nlogp(xi;θ)
Q ( θ , θ t ) − H ( θ , θ t ) Q(\theta,\theta^t)-H(\theta,\theta^t) Q(θ,θt)H(θ,θt)就是对数似然,则有
∑ i = 1 n log ⁡ p ( x i ; θ t + 1 ) − ∑ i = 1 n log ⁡ p ( x i ; θ t ) = [ Q ( θ t + 1 , θ t ) − Q ( θ t , θ t ) ] − [ H ( θ t + 1 , θ t ) − H ( θ t , θ t ) ] ( 式 1.6 ) \begin{aligned} \sum_{i=1}^n \log p(x_i;\theta^{t+1})-\sum_{i=1}^n \log p(x_i;\theta^{t})=[Q(\theta^{t+1},\theta^t)-Q(\theta^t,\theta^t)]-[H(\theta^{t+1},\theta^t)-H(\theta^t,\theta^t)](式1.6) \end{aligned} i=1nlogp(xi;θt+1)i=1nlogp(xi;θt)=[Q(θt+1,θt)Q(θt,θt)][H(θt+1,θt)H(θt,θt)](1.6)

由EM算法的M步,我们有
∑ i = 1 k ∑ j = 1 m p ( z j ∣ x i ; θ t ) l o g p ( x i , z j ; θ t + 1 ) p ( z j ∣ x i ; θ t ) ≥ ∑ i = 1 k ∑ j = 1 m p ( z j ∣ x i ; θ t ) l o g p ( x i , z j ; θ t ) p ( z j ∣ x i ; θ t ) \sum_{i=1}^k\sum_{j=1}^m p(z_j|x_i;\theta^t) log\frac{p(x_i,z_j;\theta^{t+1})}{p(z_j|x_i;\theta^t)}\geq \sum_{i=1}^k\sum_{j=1}^m p(z_j|x_i;\theta^t) log\frac{p(x_i,z_j;\theta^t)}{p(z_j|x_i;\theta^t)} i=1kj=1mp(zjxi;θt)logp(zjxi;θt)p(xi,zj;θt+1)i=1kj=1mp(zjxi;θt)logp(zjxi;θt)p(xi,zj;θt)
两端同时加上 ∑ i = 1 k ∑ j = 1 m p ( z j ∣ x i ; θ t ) log ⁡ p ( z j ∣ x i ; θ t ) \sum_{i=1}^k\sum_{j=1}^m p(z_j|x_i;\theta^t) \log p(z_j|x_i;\theta^t) i=1kj=1mp(zjxi;θt)logp(zjxi;θt)(懒,所以不写了),则有
Q ( θ t + 1 , θ t ) − Q ( θ t , θ t ) ≥ 0 Q(\theta^{t+1},\theta^t)-Q(\theta^t,\theta^t)\geq 0 Q(θt+1,θt)Q(θt,θt)0
第二项会用到Jenson不等式
H ( θ t + 1 , θ t ) − H ( θ t , θ t ) = ∑ i = 1 n ∑ j = 1 m log ⁡ p ( z i j ∣ x i ; θ t + 1 ) p ( z i j ∣ x i ; θ t ) p ( z i j ∣ x i ; θ t ) ≤ ∑ i = 1 n log ⁡ ∑ j = 1 m p ( z i j ∣ x i ; θ t + 1 ) p ( z i j ∣ x i ; θ t ) p ( z i j ∣ x i ; θ t ) = ∑ i = 1 n log ⁡ ∑ j = 1 m p ( z i j ∣ x i ; θ t + 1 ) = 0 \begin{aligned} H(\theta^{t+1},\theta^t)-H(\theta^t,\theta^t)=&\sum_{i=1}^n\sum_{j=1}^m\log \frac {p(z_{ij}|x_i;\theta^{t+1})}{p(z_{ij}|x_i;\theta^{t})}p(z_{ij}|x_i;\theta^t)\\ \leq &\sum_{i=1}^n \log \sum_{j=1}^m \frac {p(z_{ij}|x_i;\theta^{t+1})}{p(z_{ij}|x_i;\theta^{t})}p(z_{ij}|x_i;\theta^t)\\ =&\sum_{i=1}^n \log \sum_{j=1}^mp(z_{ij}|x_i;\theta^{t+1})=0 \end{aligned} H(θt+1,θt)H(θt,θt)==i=1nj=1mlogp(zijxi;θt)p(zijxi;θt+1)p(zijxi;θt)i=1nlogj=1mp(zijxi;θt)p(zijxi;θt+1)p(zijxi;θt)i=1nlogj=1mp(zijxi;θt+1)=0
其实就是E[f(x)] ≤ \leq f[E(x)],所以式1.6大于等于0,式1.5得证

至此,EM算法已经全部讲解完毕,HMM算法的推导就有用到EM算法,具体可查看机器学习——隐马尔科夫模型


EM算法应用实例

本节将采用EM算法对高斯混合模型(GMM)进行参数估计,从而利用GMM进行聚类

假定我们有一批数据,这批数据服从不同的高斯分布,这些高斯分布加权和构成了高斯混合模型(GMM),对一个数据,我们不知道它来自GMM中的哪个高斯分布,这便是隐变量,推导过程可以查看:EM算法对高斯混合模型(GMM)进行参数估计

参数估计完毕后,针对这批数据,计算隐变量的后验分布,即 p ( z ∣ x ; θ ) p(z|x;\theta) p(zx;θ),其中最大的值对应的隐变量,即为聚类结果

实验代码如下:

gmm.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# 利用均值为mu_k,方差为cov_k的正态分布计算Y的概率密度(也可以称为概率)
def phi(Y,mu_k,cov_k):
    norm=multivariate_normal(mean=mu_k,cov=cov_k)
    return norm.pdf(Y)


# EM算法中的E步
def getExpectation(Y,mu,cov,alpha):

    # 样本数,Y.shape[1]表示随机变量个数
    N=Y.shape[0]

    # 模型个数
    K=alpha.shape[0]

    # gamma用于表示样本属于各个隐变量的后验概率,即p(z|x,\theta)
    gamma=np.mat(np.zeros((N,K)))

    # 计算样本属于各个隐变量的后验概率,即p(z|x,\theta)
    prob=np.zeros((N,K))

    for k in range(K):
        prob[:,k]=phi(Y,mu[k],cov[k])
    prob=np.mat(prob)

    for k in range(K):
        gamma[:,k]=alpha[k]*prob[:,k]


    for i in range(N):
        gamma[i,:]/=np.sum(gamma[i,:])

    return gamma


# EM算法中的M步
def maximize(Y,gamma):

    # N表示样本数,D表示特征数,特征数等于一个多元高斯分布中的自变量个数
    N,D=Y.shape

    K=gamma.shape[1]

    mu=np.zeros((K,D))
    cov=[]
    alpha=np.zeros((K))

    # 更新模型中的所有参数
    for k in range(K):
        NK=np.sum(gamma[:,k])

        mu[k,:]=np.sum(np.multiply(Y,gamma[:,k]),axis=0)/NK

        cov_k=(Y-mu[k]).T*np.multiply((Y-mu[k]),gamma[:,k])/NK

        cov.append(cov_k)

        alpha[k]=NK/N

    cov=np.array(cov)
    return mu,cov,alpha


def scale_data(Y):
    for i in range(Y.shape[1]):
        max=Y[:,i].max()
        min=Y[:,i].min()
        Y[:,i]=(Y[:,i]-min)/(max-min)
    return Y

def init_params(shape,K):
    N,D=shape
    mu=np.random.rand(K,D)
    cov=np.array([np.eye(D)]*K)
    alpha=np.array([1.0/K]*K)
    return mu,cov,alpha



def GMM_EM(Y,K,times):
    Y=scale_data(Y)
    mu,cov,alpha=init_params(Y.shape,K)
    for i in range(times):
        gamma=getExpectation(Y,mu,cov,alpha)
        mu,cov,alpha=maximize(Y,gamma)
    return mu,cov,alpha

main.py

from gmm import *

Y=np.loadtxt("gmm.data")
matY=np.matrix(Y,copy=True)

K=2

mu,cov,alpha=GMM_EM(matY,K,100)

N=Y.shape[0]

gamma=getExpectation(matY,mu,cov,alpha)

category=gamma.argmax(axis=1).flatten().tolist()[0]

class1=np.array([Y[i] for i in range(N) if category[i]==0])
class2=np.array([Y[i] for i in range(N) if category[i]==1])

plt.plot(class1[:,0],class1[:,1],'rs',label="class1")
plt.plot(class2[:,0],class2[:,1],'bo',label='class2')
plt.legend(loc="best")
plt.show()

gmm.data

3.600000 79.000000
1.800000 54.000000
3.333000 74.000000
2.283000 62.000000
4.533000 85.000000
2.883000 55.000000
4.700000 88.000000
3.600000 85.000000
1.950000 51.000000
4.350000 85.000000
1.833000 54.000000
3.917000 84.000000
4.200000 78.000000
1.750000 47.000000
4.700000 83.000000
2.167000 52.000000
1.750000 62.000000
4.800000 84.000000
1.600000 52.000000
4.250000 79.000000
1.800000 51.000000
1.750000 47.000000
3.450000 78.000000
3.067000 69.000000
4.533000 74.000000
3.600000 83.000000
1.967000 55.000000
4.083000 76.000000
3.850000 78.000000
4.433000 79.000000
4.300000 73.000000
4.467000 77.000000
3.367000 66.000000
4.033000 80.000000
3.833000 74.000000
2.017000 52.000000
1.867000 48.000000
4.833000 80.000000
1.833000 59.000000
4.783000 90.000000
4.350000 80.000000
1.883000 58.000000
4.567000 84.000000
1.750000 58.000000
4.533000 73.000000
3.317000 83.000000
3.833000 64.000000
2.100000 53.000000
4.633000 82.000000
2.000000 59.000000
4.800000 75.000000
4.716000 90.000000
1.833000 54.000000
4.833000 80.000000
1.733000 54.000000
4.883000 83.000000
3.717000 71.000000
1.667000 64.000000
4.567000 77.000000
4.317000 81.000000
2.233000 59.000000
4.500000 84.000000
1.750000 48.000000
4.800000 82.000000
1.817000 60.000000
4.400000 92.000000
4.167000 78.000000
4.700000 78.000000
2.067000 65.000000
4.700000 73.000000
4.033000 82.000000
1.967000 56.000000
4.500000 79.000000
4.000000 71.000000
1.983000 62.000000
5.067000 76.000000
2.017000 60.000000
4.567000 78.000000
3.883000 76.000000
3.600000 83.000000
4.133000 75.000000
4.333000 82.000000
4.100000 70.000000
2.633000 65.000000
4.067000 73.000000
4.933000 88.000000
3.950000 76.000000
4.517000 80.000000
2.167000 48.000000
4.000000 86.000000
2.200000 60.000000
4.333000 90.000000
1.867000 50.000000
4.817000 78.000000
1.833000 63.000000
4.300000 72.000000
4.667000 84.000000
3.750000 75.000000
1.867000 51.000000
4.900000 82.000000
2.483000 62.000000
4.367000 88.000000
2.100000 49.000000
4.500000 83.000000
4.050000 81.000000
1.867000 47.000000
4.700000 84.000000
1.783000 52.000000
4.850000 86.000000
3.683000 81.000000
4.733000 75.000000
2.300000 59.000000
4.900000 89.000000
4.417000 79.000000
1.700000 59.000000
4.633000 81.000000
2.317000 50.000000
4.600000 85.000000
1.817000 59.000000
4.417000 87.000000
2.617000 53.000000
4.067000 69.000000
4.250000 77.000000
1.967000 56.000000
4.600000 88.000000
3.767000 81.000000
1.917000 45.000000
4.500000 82.000000
2.267000 55.000000
4.650000 90.000000
1.867000 45.000000
4.167000 83.000000
2.800000 56.000000
4.333000 89.000000
1.833000 46.000000
4.383000 82.000000
1.883000 51.000000
4.933000 86.000000
2.033000 53.000000
3.733000 79.000000
4.233000 81.000000
2.233000 60.000000
4.533000 82.000000
4.817000 77.000000
4.333000 76.000000
1.983000 59.000000
4.633000 80.000000
2.017000 49.000000
5.100000 96.000000
1.800000 53.000000
5.033000 77.000000
4.000000 77.000000
2.400000 65.000000
4.600000 81.000000
3.567000 71.000000
4.000000 70.000000
4.500000 81.000000
4.083000 93.000000
1.800000 53.000000
3.967000 89.000000
2.200000 45.000000
4.150000 86.000000
2.000000 58.000000
3.833000 78.000000
3.500000 66.000000
4.583000 76.000000
2.367000 63.000000
5.000000 88.000000
1.933000 52.000000
4.617000 93.000000
1.917000 49.000000
2.083000 57.000000
4.583000 77.000000
3.333000 68.000000
4.167000 81.000000
4.333000 81.000000
4.500000 73.000000
2.417000 50.000000
4.000000 85.000000
4.167000 74.000000
1.883000 55.000000
4.583000 77.000000
4.250000 83.000000
3.767000 83.000000
2.033000 51.000000
4.433000 78.000000
4.083000 84.000000
1.833000 46.000000
4.417000 83.000000
2.183000 55.000000
4.800000 81.000000
1.833000 57.000000
4.800000 76.000000
4.100000 84.000000
3.966000 77.000000
4.233000 81.000000
3.500000 87.000000
4.366000 77.000000
2.250000 51.000000
4.667000 78.000000
2.100000 60.000000
4.350000 82.000000
4.133000 91.000000
1.867000 53.000000
4.600000 78.000000
1.783000 46.000000
4.367000 77.000000
3.850000 84.000000
1.933000 49.000000
4.500000 83.000000
2.383000 71.000000
4.700000 80.000000
1.867000 49.000000
3.833000 75.000000
3.417000 64.000000
4.233000 76.000000
2.400000 53.000000
4.800000 94.000000
2.000000 55.000000
4.150000 76.000000
1.867000 50.000000
4.267000 82.000000
1.750000 54.000000
4.483000 75.000000
4.000000 78.000000
4.117000 79.000000
4.083000 78.000000
4.267000 78.000000
3.917000 70.000000
4.550000 79.000000
4.083000 70.000000
2.417000 54.000000
4.183000 86.000000
2.217000 50.000000
4.450000 90.000000
1.883000 54.000000
1.850000 54.000000
4.283000 77.000000
3.950000 79.000000
2.333000 64.000000
4.150000 75.000000
2.350000 47.000000
4.933000 86.000000
2.900000 63.000000
4.583000 85.000000
3.833000 82.000000
2.083000 57.000000
4.367000 82.000000
2.133000 67.000000
4.350000 74.000000
2.200000 54.000000
4.450000 83.000000
3.567000 73.000000
4.500000 73.000000
4.150000 88.000000
3.817000 80.000000
3.917000 71.000000
4.450000 83.000000
2.000000 56.000000
4.283000 79.000000
4.767000 78.000000
4.533000 84.000000
1.850000 58.000000
4.250000 83.000000
1.983000 43.000000
2.250000 60.000000
4.750000 75.000000
4.117000 81.000000
2.150000 46.000000
4.417000 90.000000
1.817000 46.000000
4.467000 74.000000
  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值