EM算法笔记

蓝是一种心情
88 阅读3分钟

EM算法

如果概率模型的变量都是观察变量，那么给定数据可以直接使用极大似然估计法。

当模型含有隐向量时，就不能简单地使用这些估计方法。——From 李航《统计学习方法》

如何理解上面这句话？

9.1 EM算法举例

  1. A，B，C三个硬币，出现正面的概率分别π，p和q
  2. 抛A硬币
    1. if 正面：抛B硬币 π
    2. else： 抛C硬币 (1-π)
  3. 对于抛B和C硬币的结果，如果正面记1，反面记0 （注意这里不记住A的记录）

如独立重复n=10次，观测结果：

1,1,0,1,0,0,1,0,11,1,0,1,0,0,1,0,1

现在问题：上面实验中出现1的概率是如何？

P(y∣π^,p^,q^)=P(y∣θ)P(y|\hatπ ,\hat p,\hat q)=P(y|\theta)

其中 θ\theta 为 (π,p,q)(π,p,q) 参数的估计。yy为在观测值，取值0或1。

P(y=1∣θ)=πp+(1−π)qP(y=0∣θ)=π(1−p)+(1−π)(1−q)P(y=1|\theta) = πp+(1-π)q \\ P(y=0|\theta) = π(1-p)+(1-π)(1-q)

可以将上面两个合并为一个公式：

P(y∣θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−yP(y|\theta)=πp^y(1-p)^{1-y}+(1-π)q^y(1-q)^{1-y}

现：已知观测数据Y=(Y1,Y2,...,Yn)TY=(Y_1,Y_2,..., Y_n)^T，YiY_i为独立实验的单次观测结果。出现YY序列Y1,Y2,...,YnY_1,Y_2,..., Y_n出现的概率为n个独立实验概率的积：

P(Y∣θ)=∏j=1n[P(yi∣θ)]=∏j=1n[πpy(1−p)1−y+(1−π)qy(1−q)1−y]P(Y|\theta)=\prod^n_{j=1}[P(y_i|\theta)] = \prod^n_{j=1}[πp^y(1-p)^{1-y}+(1-π)q^y(1-q)^{1-y}]

根据上面所说：如果概率模型的变量都是观察变量，那么给定数据可以直接使用极大似然估计法，我们试着给出：

比如给定数据 Y′Y' 为：

1,1,0,1,0,0,1,0,11,1,0,1,0,0,1,0,1

则出现上面序列的概率：

P(Y′∣θ)=∏j=110[P(yi∣θ)]=[πp+(1−π)q]5∗[π(1−p)+(1−π)(1−q)]5P(Y'|\theta)=\prod^{10}_{j=1}[P(y_i|\theta)] = [πp+(1-π)q]^5*[π(1-p)+(1-π)(1-q)]^5

而我们并不知道θ\theta具体取值，所以我们用极大似然估计法，取得对 θ\theta 的估计，即让上面式子取极大值时 (π,p,q)(π,p,q) 的取值 (π^,p^,q^)(\hatπ ,\hat p,\hat q)。

θ^=argmax⁡nlogP(Y∣θ)=argmax⁡[5log(πp+(1−π)q)+5log(π(1−p)+(1−π)(1−q))]=argmax⁡[5log(πp+q−πq)+5log(−πp+1+πq−q)]\hat\theta=arg\underset{n}{\max}logP(Y|\theta) \\ =arg{\max}[5log(πp+(1-π)q) + 5log(π(1-p)+(1-π)(1-q))] \\ =arg{\max}[5log(πp+q-πq) + 5log(-πp+1+πq-q)]

接下来就是求解个分量偏导=0时候的方程解，即可求出 (π^,p^,q^)(\hatπ ,\hat p,\hat q)的值。

获取三个方程：

p−qπp+q−πq=p−q−πp+1+πq−qππp+q−πq=π−πp+1+πq−q1−ππp+q−πq=1−π−πp+1+πq−q\frac{p-q}{\pi p + q - \pi q} = \frac{p - q}{- \pi p + 1 + \pi q - q} \\ \frac{\pi}{\pi p + q - \pi q} = \frac{\pi}{- \pi p + 1 + \pi q - q} \\ \frac{1 - \pi}{\pi p + q - \pi q} = \frac{1 - \pi}{- \pi p + 1 + \pi q - q}

分析可得：1. 分母相等，p/q/π互相依赖，解不出来参数。

查看链接：[5分钟学算法] #06 EM算法 你到底是哪个班级的

EM算法用于求解这种问题的迭代算法。使用到Jensen不等式：zhuanlan.zhihu.com/p/28249050/

我们要通过EM算法极大化：

链式法则：链式法则用于计算多个事件的联合概率，可以表示为 P(A1,...,An)=P(A1)P(A2∣A1)P(A3∣A1A2)...P(An∣A1A2...An−1)P(A1,...,An)=P(A1)P(A2∣A1)P(A3∣A1A2)...P(An∣A1A2...A{n−1})

全概率公式：如果事件B1, B2, ..., Bn构成一个完备事件组（即它们是互斥的，并且它们的并集是整个样本空间），那么对于任意事件A，有

P(A)=∑i=1nP(A∣Bi)P(Bi)=∑i=1nP(A,Bi)P(A)=∑_{i=1}^nP(A∣Bi)P(Bi) = ∑_{i=1}^nP(A,B_i)

利用上面法则以及贝叶斯公式可以得到：

L(θ)=logP(Y∣θ)=log∑ZP(Y,Z∣θ)=log(∑ZP(Y∣Z,θ))P(Z∣θ)L(\theta) = logP(Y|\theta) = log\sum_ZP(Y,Z|\theta) \\ =log(\sum_ZP(Y|Z,\theta))P(Z|\theta)

EM会对θ\theta进行更新，希望L(θ)L(\theta)能不断增大直至收敛得到极大值。即L(θ)−L(θ(i))L(\theta)-L(\theta^{(i)})，随着i的增大，两者差越来越小。

L(θ)−L(θ(i))=logP(Y∣θ)−logP(Y∣θ(i))=log(∑ZP(Y∣Z,θ))P(Z∣θ)−logP(Y∣θ(i))=log(∑ZP(Z∣Y,θ(i))P(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(i)))−logP(Y∣θ(i))L(\theta)-L(\theta^{(i)}) = logP(Y|\theta) - logP(Y|\theta^{(i)})\\ =log(\sum_ZP(Y|Z,\theta))P(Z|\theta) - logP(Y|\theta^{(i)}) \\ =log(\sum_ZP(Z|Y,\theta^{(i)}) \frac{P(Y|Z,\theta))P(Z|\theta)}{P(Z|Y,\theta^{(i)})})-logP(Y|\theta^{(i)})

根据Jensen不等式，对于函数为凹函数的 loglog 存在：

log(∑jλjyj)≥∑jλjlog(yi)，其中λj≥0,∑jλj=1log(\sum_j\lambda_jy_j)\ge\sum_j\lambda_jlog(y_i) \hspace{15pt} ，其中 \lambda_j\ge0,\sum_j\lambda_j=1

上式子继续化简：

≥∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(1))−logP(Y∣θ(i))=∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(i))P(Y∣θ(i))\ge\sum_ZP(Z|Y,\theta^{(i)})log \frac{P(Y|Z,\theta))P(Z|\theta)}{P(Z|Y,\theta^{(1)})}-logP(Y|\theta^{(i)})\\ =\sum_ZP(Z|Y,\theta^{(i)})log \frac{P(Y|Z,\theta))P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}

所以有：

L(θ)−L(θ(i))≥∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(i))P(Y∣θ(i))L(\theta)-L(\theta^{(i)}) \ge\sum_ZP(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))+MAX(∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(i))P(Y∣θ(i)))L(\theta^{(i+1)}) = L(\theta^{(i)}) +MAX(\sum_ZP(Z|Y,\theta^{(i)})log \frac{P(Y|Z,\theta))P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})})

任务则去求：

argmax⁡θ(∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)P(Z∣Y,θ(i))P(Y∣θ(i)))arg\underset{\theta}{\max}(\sum_ZP(Z|Y,\theta^{(i)})log \frac{P(Y|Z,\theta))P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})})

去掉对θ\theta而言是常数的项目，得：

θ(i+1)=argmax⁡θ(∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)=argmax⁡θQ(θ,θ(i))\theta^{(i+1)} =arg\underset{\theta}{\max}(\sum_ZP(Z|Y,\theta^{(i)})log {P(Y|Z,\theta))P(Z|\theta)} \\ =arg\underset{\theta}{\max} Q(\theta, \theta^{(i)})

E步（求期望）：

即：似然函数log(P(Y,Z∣θ))log(P(Y,Z|\theta))关于给定观测数据Y和当前参数θ(i)\theta^{(i)}对未观测数据Z的条件概率分布的期望P(Z∣Y,θ(i))P(Z|Y,\theta^{(i)})，记为Q(θ,θ(i))Q(\theta, \theta^{(i)}).

Q(θ,θ(i))=∑ZP(Z∣Y,θ(i))logP(Y∣Z,θ))P(Z∣θ)=∑ZP(Z∣Y,θ(i))logP(Y,Z∣θ)Q(\theta, \theta^{(i)}) = \sum_ZP(Z|Y,\theta^{(i)})log {P(Y|Z,\theta))P(Z|\theta)}\\ = \sum_ZP(Z|Y,\theta^{(i)})log {P(Y,Z|\theta)}

我们看第一部分，利用贝叶斯公式：

1）对于选了硬币B

P(zj=B∣yj,θ(i))=P(yj∣zj,θ(i))P(zj∣θ(i))P(yj∣θ(i))=π(i)[p(i)]yj(1−p(i))1−yjπ(i)[p(i)]yj(1−p(i))1−yj+(1−π(i))[q(i)]yj(1−q(i))1−yj≜ujP(z_j=B|y_j,\theta^{(i)}) =\frac{P(y_j|z_j,\theta^{(i)}) P(z_j|\theta^{(i)})}{P(y_j|\theta^{(i)})} \\ =\frac{π^{(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}} \triangleq u_j

2）对于选了硬币C

P(zj=C∣yj,θ(i))=P(yj∣zj,θ(i))P(zj∣θ(i))P(yj∣θ(i))=(1−π(i))[q(i)]yj(1−q(i))1−yjπ(i)[p(i)]yj(1−p(i))1−yj+(1−π(i))[q(i)]yj(1−q(i))1−yj=1−ujP(z_j=C|y_j,\theta^{(i)}) =\frac{P(y_j|z_j,\theta^{(i)}) P(z_j|\theta^{(i)})}{P(y_j|\theta^{(i)})} \\ =\frac{(1-π^{(i)})[q^{(i)}]^{y_j}(1-q^{(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}}=1- u_j

那么：

Q(θ,θ(i))=∑ZP(Z∣Y,θ(i))logP(Y,Z∣θ)=∑j=1n(P(zj=B∣yj,θ(i))logP(yj,zj=B∣θ)+P(zj=C∣yj,θ(i))logP(yj,zj=C∣θ))=∑j=1n(ujlog(πpjyjpj1−yj)+(1−uj)log((1−π)qjyjqj(1−yj)))Q(\theta, \theta^{(i)}) = \sum_ZP(Z|Y,\theta^{(i)})log {P(Y,Z|\theta)} \\ =\sum_{j=1}^{n}(P(z_j=B|y_j,\theta^{(i)})log {P(y_j,z_j=B|\theta)} +P(z_j=C|y_j,\theta^{(i)}) log {P(y_j,z_j=C|\theta)}) \\ =\sum_{j=1}^{n}(u_jlog(πp^{y_j}_{j}p_j^{1-y_j})+(1-u_j)log((1-π)q_j^{y_j}q_j^{(1-y_j)}))

M步（求极大化）：

1） 对π求偏导：

因为uju_j中都是预估值，所以跟实际参数无关，可以作为常数处理。

∂Q∂π=∑j=1n(ujpjyjpj1−yjπpjyjpj1−yj+(1−uj)−qjyjqj(1−yj)(1−π)qjyjqj(1−yj))=∑j=1n(ujπ−(1−uj)1−π)=∑j=1nuj−π(1−π)π\frac{\partial Q}{\partial π} = \sum_{j=1}^{n}(u_j \frac{p^{y_j}_{j}p_j^{1-y_j}}{πp^{y_j}_{j}p_j^{1-y_j}}+(1-u_j)\frac{-q_j^{y_j}q_j^{(1-y_j)}}{(1-π)q_j^{y_j}q_j^{(1-y_j)}})\\ =\sum_{j=1}^{n}(\frac{u_j}{π}-\frac{(1-u_j)}{1-π})\\ =\sum_{j=1}^{n}\frac{u_j-π}{(1-π)π}

令上式等于0，因为π!=0π!=0，所以：

∂Q∂π=∑j=1n(uj−π)=∑j=1nuj−∑j=1π=∑j=1nuj−nπ=0\frac{\partial Q}{\partial π} =\sum_{j=1}^{n}{(u_j-π)}\\ =\sum_{j=1}^{n}u_j-\sum_{j=1}π\\ =\sum_{j=1}^{n}u_j -nπ=0

则，

π=∑j=1nujnπ =\frac{\sum_{j=1}^{n}u_j}{n}

作为下一次迭代的估计参数。

2）对p求偏导：

p=∑j=1nujyj∑j=1nujp = \frac{\sum_{j=1}^nu_jy_j}{\sum_{j=1}^nu_j}

3）对q求偏导：

q=∑j=1n(1−uj)yj∑j=1n(1−uj)q =\frac{\sum_{j=1}^n(1-u_j)y_j}{\sum_{j=1}^n(1-u_j)}

总结：

E：求的是Q函数，他是一个期望

M：对Q函数采用极大似然方法估计参数

avatar
文章
阅读
粉丝