算法|计算建模之EM算法机器学习|概率论

EM(Expectation-Maximization)算法于1977年提出，被用于解决如下的参数估计问题：
已知（或者至少我们的模型如此假定）随机变量 X X X服从概率密度函数 p ( X , Z ; θ ) , p(X,Z; \theta), p(X,Z; θ),其中 θ \theta θ为待确定的参数， Z Z Z为隐变量(latent variable)，是由于某些原因我们无法直接观测得到的随机变量。EM算法可以在我们只能观测到 X X X的情况下对参数 θ \theta θ进行估计。下面我们举两个例子来说明该算法。
例一这个经典的例子出自1,问题描述如下：
我们有两枚硬币A和B，它们投出正面的概率分别为 θ A \theta_A θA?和 θ B , \theta_B, θB?,是待估计的参数。假如我们现在投掷 5 5 5轮，每轮先以等概率选择其中一枚，再投掷这枚硬币 10 10 10次，得到如下结果：

轮次	硬币	结果
1	B	5正5反
2	A	9正1反
3	A	8正2反
4	B	4正6反
5	A	7正3反

这个模型中，每轮的结果（即观测变量）可以正面出现的次数唯一确定，记为 y y y; 此外，对结果有影响的变量还有每轮选择的硬币类型(A/B)，我们记为 z z z,当 z = 1 z=1 z=1时，表示我们选择了 A A A硬币；否则当 z = 0 z=0 z=0时，表示此时选择 B B B硬币。
如果这时让我们估算未知的参数 θ A \theta_A θA?，可以直接用轮次2,3,5的正面次数除以总次数( θ A ^ = 9 + 8 + 7 30 = 4 5 \hat{\theta_A}=\frac{9+8+7}{30}=\frac45 θA?^?=309+8+7?=54?); 但是为了与后面的EM算法的思想保持一致，我们稍微绕一下远：可以写出如下的概率分布函数
P ( y , z ; θ ) = P ( y ∣ z ; θ ) P ( z ) = 1 2 P ( y ∣ z ) = { 1 2 C 10 y θ A y ( 1 ? θ A ) 10 ? y , z = 1 1 2 C 10 y θ B y ( 1 ? θ B ) 10 ? y , z = 0 P(y,z; \pmb\theta)=P(y|z; \pmb\theta)P(z)=\frac12P(y|z)=\left\{ \begin{aligned} \frac12C_{10}^{y}\theta_A^y(1-\theta_A)^{10-y}, z=1\\ \frac12C_{10}^{y}\theta_B^y(1-\theta_B)^{10-y}, z=0 \end{aligned} \right. P(y,z; θθ)=P(y∣z; θθ)P(z)=21?P(y∣z)=? ? ??21?C10y?θAy?(1?θA?)10?y,z=121?C10y?θBy?(1?θB?)10?y,z=0?
容易知道它等价于下面的式子：
P ( y , z ; θ ) = 1 2 z C 10 y θ A y ( 1 ? θ A ) 10 ? y + 1 2 ( 1 ? z ) C 10 y θ B y ( 1 ? θ B ) 10 ? y P(y,z; \pmb \theta)=\frac12zC_{10}^{y}\theta_A^y(1-\theta_A)^{10-y}+\frac12(1-z)C_{10}^{y}\theta_B^y(1-\theta_B)^{10-y} P(y,z; θθ)=21?zC10y?θAy?(1?θA?)10?y+21?(1?z)C10y?θBy?(1?θB?)10?y
其中 z ∈ { 0 , 1 } , y ∈ { 0 , 1 , . . . , n } , θ = ( θ A , θ B ) . z\in\{0,1\},y\in\{0,1,...,n\},\pmb\theta=(\theta_A,\theta_B). z∈{0,1},y∈{0,1,...,n},θθ=(θA?,θB?).
为了估算式子中的 θ \pmb\theta θθ的值，可以用概率论中的极大似然估计构造对数似然函数：
L L ( θ ∣ y , z ) = l n ∏ i = 1 5 P ( y i , z i ; θ ) = l n ∏ i = 1 5 [ 1 2 z i C 10 y i θ A y i ( 1 ? θ A ) 10 ? y i + 1 2 ( 1 ? z i ) C 10 y i θ B y i ( 1 ? θ B ) 10 ? y i ] = ∑ i = 1 5 l n [ 1 2 z i C 10 y i θ A y i ( 1 ? θ A ) 10 ? y i + 1 2 ( 1 ? z i ) C 10 y i θ B y i ( 1 ? θ B ) 10 ? y i ] \begin{aligned}LL(\pmb\theta|y,z) & =ln\prod\limits_{i=1}^5P(y_i,z_i; \pmb\theta)\\ & =ln\prod\limits_{i=1}^{5}[\frac12z_iC_{10}^{y_i}\theta_A^{y_i}(1-\theta_A)^{10-y_i}+\frac12(1-z_i)C_{10}^{y_i}\theta_B^{y_i}(1-\theta_B)^{10-y_i}]\\ & =\sum\limits_{i=1}^{5}ln[\frac12z_iC_{10}^{y_i}\theta_A^{y_i}(1-\theta_A)^{10-y_i}+\frac12(1-z_i)C_{10}^{y_i}\theta_B^{y_i}(1-\theta_B)^{10-y_i}] \\ \end{aligned} LL(θθ∣y,z)?=lni=1∏5?P(yi?,zi?; θθ)=lni=1∏5?[21?zi?C10yi??θAyi??(1?θA?)10?yi?+21?(1?zi?)C10yi??θByi??(1?θB?)10?yi?]=i=1∑5?ln[21?zi?C10yi??θAyi??(1?θA?)10?yi?+21?(1?zi?)C10yi??θByi??(1?θB?)10?yi?]?
（注：对数似然函数 L L ( θ ∣ y , z ) LL(\pmb\theta|y,z) LL(θθ∣y,z)的记号表示 y , z y,z y,z是已知的，此时 θ \pmb\theta θθ才是该函数的自变量）
之后可以令该函数对 θ \pmb\theta θθ偏导为0，以求出 θ \pmb\theta θθ的值使得 L L ( θ ∣ y , z ) LL(\pmb\theta|y,z) LL(θθ∣y,z)最大，即令 ? ( L L ) / ? θ = 0 , \partial (LL)/ \partial{\pmb\theta}=0, ?(LL)/?θθ=0,解出 θ \pmb\theta θθ的值即可。
这个问题的难度可以升一下级：当我们不知道每轮取到的硬币时（见下表），能否估计参数 θ \pmb\theta θθ的值呢？EM算法就可以用来解决这个问题。

轮次	硬币	结果
1	？	5正5反
2	？	9正1反
3	？	8正2反
4	？	4正6反
5	？	7正3反

不知道每轮取到的硬币，也就是上面的概率分布函数 P ( y , z ) P(y,z) P(y,z)中， z z z的值无法确定了，从而在求解对数似然函数偏导数为 0 0 0时，我们没法解出 θ \pmb\theta θθ的确定值（因为 z i z_i zi?是未知的）。EM算法的解决方案是，初始化一个参数 θ 0 , \pmb\theta^0, θθ0,于是对于每个样本 y i y_i yi?,概率分布函数
P ( z ∣ y , θ 0 ) = 1 2 z C 10 y θ A 0y ( 1 ? θ A 0 ) 10 ? y + 1 2 ( 1 ? z ) C 10 y θ B 0y ( 1 ? θ B 0 ) 10 ? y P(z|y,\pmb \theta^0)=\frac12zC_{10}^{y}\theta_A^{0\ y}(1-\theta_A^0)^{10-y}+\frac12(1-z)C_{10}^{y}\theta_B^{0\ y}(1-\theta^0_B)^{10-y} P(z∣y,θθ0)=21?zC10y?θA0 y?(1?θA0?)10?y+21?(1?z)C10y?θB0 y?(1?θB0?)10?y
中，只有 z i z_i zi?是不确定的了（这里的符号有一些变化，即观测变量 y i y_i yi?变成了条件，已知），也就是说，我们可以求出 z i z_i zi?的概率分布。那么，我们就可以求出对数似然函数 L L ( θ ∣ y , z ) LL(\pmb\theta|y,z) LL(θθ∣y,z)的期望 E z [ L L ( θ ∣ y , z ) ] \mathbb{E}_z[LL(\pmb\theta|y,z)] Ez?[LL(θθ∣y,z)]，来代替原来的对数似然函数（取期望这个操作消去了不确定的 z z z，使其变为仅与 θ \pmb\theta θθ有关的函数）。这时我们再求出使得 E z [ L L ( θ ∣ y , z ) ] \mathbb{E}_z[LL(\pmb\theta|y,z)] Ez?[LL(θθ∣y,z)]最大的 θ 1 \pmb\theta^1 θθ1，作为新一轮的 θ \pmb\theta θθ的估计值，重复上面的操作直到参数 θ \pmb\theta θθ收敛（收敛性在此不证明，分析可参见2），算法停止。
总结一下EM算法的步骤：
(0)：用随机数初始化位置参数 θ 0 \pmb\theta^0 θθ0;
(1)（E步）将观测变量 y i y_i yi?代入 L L ( θ ∣ y , z ) = l n ∏ i = 1 N P ( y i , z i ; θ ) LL(\pmb\theta|y,z)=ln\prod\limits_{i=1}^NP(y_i,z_i; \pmb \theta) LL(θθ∣y,z)=lni=1∏N?P(yi?,zi?; θθ),求出对数似然函数关于 z z z的期望 E z [ L L ( θ ∣ y , z ) ] \mathbb{E}_z[LL(\pmb\theta|y,z)] Ez?[LL(θθ∣y,z)];
(2)（M步）求出使得(1)中期望最大的参数 θ t + 1 . \pmb\theta^{t+1}. θθt+1.这里， z z z服从的概率分布由上一轮的参数 θ t \pmb\theta^t θθt确定。
(3) 若迭代次数到达指定上限或 θ \pmb\theta θθ已经收敛（即 θ t + 1 \pmb\theta^{t+1} θθt+1和 θ t \pmb\theta^t θθt的差距已经非常小），结束算法；否则回到(1)继续进行迭代。
说明：这里出现的 θ \pmb\theta θθ和 θ t \pmb\theta^t θθt容易让人弄混。我们在EM算法中的目标函数是 E z [ L L ( θ ∣ y , z ) ] \mathbb{E}_z[LL(\pmb\theta|y,z)] Ez?[LL(θθ∣y,z)]，这是一个关于未知参数 θ \pmb\theta θθ的函数；估计出来的 θ t \pmb\theta^{t} θθt用于确定隐变量 z z z的分布，以计算 L L ( θ ∣ y , z ) LL(\pmb\theta|y,z) LL(θθ∣y,z)关于 z z z的期望。
下面我们对上述硬币问题应用EM算法，来熟悉一下这个过程。我们这里把轮次记为 N N N，每轮投掷次数为 n n n（上面的例子中， N = 5 , n = 10 N=5,n=10 N=5,n=10）.设每轮观测到正面朝上的次数为 y i , y_i, yi?,每轮选出的硬币为 z i ∈ { 0 , 1 } z_i\in\{0,1\} zi?∈{0,1},是一个隐变量。
（1）初始化 θ 0 , \pmb\theta^0, θθ0,可以调用随机函数生成，这里不指定具体数值。
（2）写出概率分布函数：
P ( y , z ) = 1 2 z C n y θ A y ( 1 ? θ A ) n ? y + 1 2 ( 1 ? z ) C n y θ B y ( 1 ? θ B ) n ? y P(y,z)=\frac12zC_{n}^{y}\theta_A^y(1-\theta_A)^{n-y}+\frac12(1-z)C_{n}^{y}\theta_B^y(1-\theta_B)^{n-y} P(y,z)=21?zCny?θAy?(1?θA?)n?y+21?(1?z)Cny?θBy?(1?θB?)n?y
（3）表示对数似然函数关于隐变量 z z z的期望：
E z [ L L ( θ ∣ y , z ) ] = E z { ∑ i = 1 N l n [ 1 2 z i C n y i θ A y i ( 1 ? θ A ) n ? y i + 1 2 ( 1 ? z i ) C n y i θ B y ( 1 ? θ B ) n ? y i ] } = ∑ i = 1 N E z { l n [ 1 2 z i C n y i θ A y i ( 1 ? θ A ) n ? y i + 1 2 ( 1 ? z i ) C n y i θ B y ( 1 ? θ B ) n ? y i ] } ( 把期望号放进求和号 ) = ∑ i = 1 N { P ( z = 1 ∣ y i , θ t ) l n [ 1 2 C n y i θ A y i ( 1 ? θ A ) n ? y i ] + P ( z = 0 ∣ y i , θ t ) l n [ C n y i θ B y ( 1 ? θ B ) n ? y i ] } ( 离散型随机变量期望的定义 ) \begin{aligned}\mathbb{E_z}[LL(\pmb\theta|y,z)] & = \mathbb{E_z}\{\sum\limits_{i=1}^Nln[\frac12z_iC_{n}^{y_i}\theta_A^{y_i}(1-\theta_A)^{n-y_i}+\frac12(1-z_i)C_{n}^{y_i}\theta_B^y(1-\theta_B)^{n-y_i}]\} \\ & = \sum\limits_{i=1}^N\mathbb{E_z}\{ln[\frac12z_iC_{n}^{y_i}\theta_A^{y_i}(1-\theta_A)^{n-y_i}+\frac12(1-z_i)C_{n}^{y_i}\theta_B^y(1-\theta_B)^{n-y_i}]\} (把期望号放进求和号)\\ & =\sum\limits_{i=1}^N\{P(z=1|y_i,\pmb\theta^t)ln[\frac12C_{n}^{y_i}\theta_A^{y_i}(1-\theta_A)^{n-y_i}]+P(z=0|y_i,\pmb\theta^t)ln[C_{n}^{y_i}\theta_B^y(1-\theta_B)^{n-y_i}] \}(离散型随机变量期望的定义)\\ \end{aligned} Ez?[LL(θθ∣y,z)]?=Ez?{i=1∑N?ln[21?zi?Cnyi??θAyi??(1?θA?)n?yi?+21?(1?zi?)Cnyi??θBy?(1?θB?)n?yi?]}=i=1∑N?Ez?{ln[21?zi?Cnyi??θAyi??(1?θA?)n?yi?+21?(1?zi?)Cnyi??θBy?(1?θB?)n?yi?]}(把期望号放进求和号)=i=1∑N?{P(z=1∣yi?,θθt)ln[21?Cnyi??θAyi??(1?θA?)n?yi?]+P(z=0∣yi?,θθt)ln[Cnyi??θBy?(1?θB?)n?yi?]}(离散型随机变量期望的定义)?
（4）求使得 E z [ L L ( θ ∣ y , z ) ] \mathbb{E_z}[LL(\pmb\theta|y,z)] Ez?[LL(θθ∣y,z)]最大的 θ \pmb\theta θθ，需要对 θ A , θ B \theta_A,\theta_B θA?,θB?求偏导,这里以 θ A \theta_A θA?为例：
? E z ? θ A = ∑ i = 1 N P ( z = 1 ∣ y i , θ t ) 1 2 C n y i θ A y i ? 1 ( 1 ? θ A ) n ? y i ? 1 ( y i ? n θ A ) 1 2 C n y i θ A y i ( 1 ? θ A ) n ? y i = ∑ i = 1 N P ( z = 1 ∣ y i , θ t ) y i ? n θ A θ A ( 1 ? θ A ) \begin{aligned}\frac{\partial\mathbb{E_z}}{\partial\theta_A} & = \sum\limits_{i=1}^NP(z=1|y_i,\pmb\theta^t)\frac{\frac12C_n^{y_i}\theta_A^{y_i-1}(1-\theta_A)^{n-y_i-1}(y_i-n\theta_A)}{\frac12C_{n}^{y_i}\theta_A^{y_i}(1-\theta_A)^{n-y_i}}\\ & = \sum\limits_{i=1}^NP(z=1|y_i,\pmb\theta^t)\frac{y_i-n\theta_A}{\theta_A(1-\theta_A)}\end{aligned} ?θA??Ez???=i=1∑N?P(z=1∣yi?,θθt)21?Cnyi??θAyi??(1?θA?)n?yi?21?Cnyi??θAyi??1?(1?θA?)n?yi??1(yi??nθA?)?=i=1∑N?P(z=1∣yi?,θθt)θA?(1?θA?)yi??nθA???
令 ? E z ? θ A = 0 , \frac{\partial\mathbb{E_z}}{\partial\theta_A}=0, ?θA??Ez??=0,
有
∑ i = 1 N P ( z = 1 ∣ y i , θ t ) ( y i ? n θ A ) = 0 , ∑ i = 1 N n θ A P ( z = 1 ∣ y i , θ t ) = ∑ i = 1 N y i P ( z = 1 ∣ y i , θ t ) \sum\limits_{i=1}^NP(z=1|y_i,\pmb\theta^t)(y_i-n\theta_A)=0,\\ \sum\limits_{i=1}^Nn\theta_AP(z=1|y_i,\pmb\theta^t)=\sum\limits_{i=1}^Ny_iP(z=1|y_i,\pmb\theta^t) i=1∑N?P(z=1∣yi?,θθt)(yi??nθA?)=0,i=1∑N?nθA?P(z=1∣yi?,θθt)=i=1∑N?yi?P(z=1∣yi?,θθt)
解得
θ A t + 1 = 1 n ? ∑ i = 1 N y i ? P ( z = 1 ∣ y i , θ t ) ∑ i = 1 N P ( z = 1 ∣ y i , θ t )( 1 ) \theta_A^{t+1}=\frac1n\cdot \frac{\sum\limits_{i=1}^Ny_i\cdot P(z=1|y_i,\pmb\theta^t)}{\sum\limits_{i=1}^NP(z=1|y_i,\pmb\theta^t)}\ \ \ \ \ \ \ \ \ \ \ (1) θAt+1?=n1??i=1∑N?P(z=1∣yi?,θθt)i=1∑N?yi??P(z=1∣yi?,θθt)?(1)
同理
θ B t + 1 = 1 n ? ∑ i = 1 N y i ? P ( z = 0 ∣ y i , θ t ) ∑ i = 1 N P ( z = 0 ∣ y i , θ t )( 2 ) \theta_B^{t+1}=\frac1n\cdot \frac{\sum\limits_{i=1}^Ny_i\cdot P(z=0|y_i,\pmb\theta^t)}{\sum\limits_{i=1}^NP(z=0|y_i,\pmb\theta^t)}\ \ \ \ \ \ \ \ \ \ \ (2) θBt+1?=n1??i=1∑N?P(z=0∣yi?,θθt)i=1∑N?yi??P(z=0∣yi?,θθt)?(2)
其中(注意 θ A , θ B \theta_A,\theta_B θA?,θB?上面的角标 t t t)
P ( z = 1 ∣ y i , θ t ) = C n y i θ A ty i ( 1 ? θ A t ) n ? y i C n y i θ A ty i ( 1 ? θ A t ) n ? y i + C n y i θ B ty i ( 1 ? θ B t ) n ? y i ( 已知 y i 次正面朝上，投的硬币是 A 的概率 ) = θ A ty i ( 1 ? θ A t ) n ? y i θ A ty i ( 1 ? θ A t ) n ? y i + θ B ty i ( 1 ? θ B t ) n ? y i . \begin{aligned}P(z=1|y_i,\pmb\theta^t) & =\frac{C_n^{y_i}\theta_A^{t\ y_i}(1-\theta_A^t)^{n-y_i}}{C_n^{y_i}\theta_A^{t\ y_i}(1-\theta_A^t)^{n-y_i}+C_n^{y_i}\theta_B^{t\ y_i}(1-\theta_B^t)^{n-y_i}}(已知y_i次正面朝上，投的硬币是A的概率) \\ & =\frac{\theta_A^{t\ y_i}(1-\theta_A^t)^{n-y_i}}{\theta_A^{t\ y_i}(1-\theta_A^t)^{n-y_i}+\theta_B^{t\ y_i}(1-\theta_B^t)^{n-y_i}}. \end{aligned} P(z=1∣yi?,θθt)?=Cnyi??θAt yi??(1?θAt?)n?yi?+Cnyi??θBt yi??(1?θBt?)n?yi?Cnyi??θAt yi??(1?θAt?)n?yi??(已知yi?次正面朝上，投的硬币是A的概率)=θAt yi??(1?θAt?)n?yi?+θBt yi??(1?θBt?)n?yi?θAt yi??(1?θAt?)n?yi??.?
那么
P ( z = 0 ∣ y i , θ t ) = 1 ? P ( z = 1 ∣ y i , θ t ) ( 对立事件 ) = θ B ty i ( 1 ? θ B t ) n ? y i θ A ty i ( 1 ? θ A t ) n ? y i + θ B ty i ( 1 ? θ B t ) n ? y i . \begin{aligned}P(z=0|y_i,\pmb\theta^t) & =1-P(z=1|y_i,\pmb\theta^t)(对立事件)\\ & = \frac{\theta_B^{t\ y_i}(1-\theta_B^t)^{n-y_i}}{\theta_A^{t\ y_i}(1-\theta_A^t)^{n-y_i}+\theta_B^{t\ y_i}(1-\theta_B^t)^{n-y_i}}.\end{aligned} P(z=0∣yi?,θθt)?=1?P(z=1∣yi?,θθt)(对立事件)=θAt yi??(1?θAt?)n?yi?+θBt yi??(1?θBt?)n?yi?θBt yi??(1?θBt?)n?yi??.?
因此我们可以用式 ( 1 ) ( 2 ) (1)(2) (1)(2)迭代地更新参数直到收敛。下面我们用python模拟一下这个过程，验证式子的正确性以及算法的收敛性：

import random import numpy as np# 每轮的投掷次数n n = 10 thetaA, thetaB = 0.5, 0.7# 模拟概率为p，次数为n的二项分布(投硬币)，返回正面朝上的次数 def observe(n, p): sum = 0 for i in range(n): t = random.random() if t < p: sum += 1 return sum# N 为数据集大小，返回 N 轮后的结果 def get_dataset(N): dataset = [] for i in range(N): if random.random() < 0.5: dataset.append(observe(n, thetaA)) else: dataset.append(observe(n, thetaB)) return dataset# 对dataset应用EM算法 def EM(dataset) # y是N维向量，表示观测结果之集 y = np.array(dataset, dtype = 'float') # 初始化估计参数。random.random()返回0~1中的小数 thetaA_hat, thetaB_hat = random.random(), random.random()iteration = 0 while True: iteration += 1 A = ((thetaA_hat ** y) * (1 - thetaA_hat) ** (n - y)) B = ((thetaB_hat ** y) * (1 - thetaB_hat) ** (n - y)) # f,g 分别表示上文的P(z=1|y,theta),P(z=0|y,theta),是N维向量 # 向量的第 i 维代表P(z|y_i,theta) f = A / (A + B) g = B / (A + B)# 这个更新与上文(1)(2)本质相同，不过用了numpy中的点乘简化运算# 即 [y_1, y_2, ..., y_N] 与 [P(z=1|y_1,theta), ..., P(z=1|y_N,theta)]的点乘 new_thetaA_hat = np.dot(y, f) / n / np.sum(f)# 即 [y_1, y_2, ..., y_N] 与 [P(z=0|y_1,theta), ..., P(z=0|y_N,theta)]的点乘 new_thetaB_hat = np.dot(y, g) / n / np.sum(g)# 计算一下两次参数之间的变化量(欧几里得距离平方) conv = (thetaA_hat - new_thetaA_hat) ** 2 + (thetaB_hat - new_thetaB_hat) ** 2 if conv < 1e-10: breakthetaA_hat = new_thetaA_hat thetaB_hat = new_thetaB_hat print(f'iteration = {iteration}') return thetaA_hat, thetaB_hat # 做10000轮实验 dataset = get_dataset(10000) # 重复调用20次EM算法，查看每次是否收敛到一个结果，以及迭代次数 for i in range(20): print(EM(dataset))

运行结果：

iteration = 15 (0.4977663316100947, 0.6962688553151227) iteration = 15 (0.4977649567622754, 0.6962702416641224) iteration = 16 (0.497762899611069, 0.6962724001030222) iteration = 15 (0.6962733823360499, 0.4977619579412015) iteration = 14 (0.49776333638894465, 0.6962718926393595) # ...(后面省略)

可以看到概率值与真实值 ( 0.5 , 0.7 ) (0.5,0.7) (0.5,0.7)非常接近了，平均只需要迭代十几次即可收敛。但我们发现有的时候结果会收敛于 ( 0.7 , 0.5 ) (0.7,0.5) (0.7,0.5)，与我们初始设置的参数相反。这其实可以理解：我们已经“屏蔽”了每次投硬币的编号，所以我们至多只能得出“有两枚硬币，其中一枚有0.5的几率朝上，另一枚有0.7的几率朝上”的结论。至于它们的编号是无关紧要也无从得知的。
例二这个例子出自3,也是哈工大课件中的例子，问题描述如下：
在一个模式识别任务中，物体大致上可分为两类：暗色物体和亮色物体。其中暗色物体还可以分为方形和圆形。一个特征提取器可以将物体分为暗色物体 ( y 1 ) (y_1) (y1?)和亮色物体 ( y 2 ) (y_2) (y2?)，却无法区分方形和圆形的物体。此外，我们已知这三种物体（暗圆 x 1 x_1 x1?，暗方 x 2 x_2 x2?，亮色 x 3 x_3 x3?）服从如下的三项分布，它是二项分布的推广：
P ( x 1 , x 2 , x 3 ; θ ) = ( x 1 + x 2 + x 3 ) ! x 1 ! x 2 ! x 3 ! ( 1 4 ) x 1 ( 1 + θ 4 ) x 2 ( 2 ? θ 4 ) x 3θ ∈ ( ? 1 , 2 ) . P(x_1,x_2,x_3; \theta)=\frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}(\frac14)^{x_1}(\frac{1+\theta}{4})^{x_2}(\frac{2-\theta}{4})^{x_3}\ \ \ \ \ \ \theta \in (-1,2). P(x1?,x2?,x3?; θ)=x1?!x2?!x3?!(x1?+x2?+x3?)!?(41?)x1?(41+θ?)x2?(42?θ?)x3?θ∈(?1,2).
在这个问题中，观测变量为 Y = ( y 1 , y 2 ) , Y=(y_1,y_2), Y=(y1?,y2?),全部变量为 X = ( x 1 , x 2 , x 3 ) X=(x_1,x_2,x_3) X=(x1?,x2?,x3?)，有
y 1 = x 1 + x 2 y 2 = x 3 y_1=x_1+x_2\\ y_2=x_3 y1?=x1?+x2?y2?=x3?
隐变量为 x 1 , x 2 x_1,x_2 x1?,x2?（ x 3 = y 2 x_3=y_2 x3?=y2?，因此是可以被观测的）。我们尝试对这个问题应用EM算法，以估计概率密度函数中的未知参数 θ \theta θ。设 X i = ( x i 1 , x i 2 , x i 3 ) X_i=(x_{i1},x_{i2},x_{i3}) Xi?=(xi1?,xi2?,xi3?)为第 i i i个样本的全部变量， Y i = ( y i 1 , y i 2 ) Y_i=(y_{i1},y_{i2}) Yi?=(yi1?,yi2?)为第 i i i个样本的观测变量，EM算法流程如下：
（1）初始化 θ 0 \pmb\theta^0 θθ0;
（2）表示对数似然函数关于隐变量 X h X_h Xh?的期望(由于我们可以求出 x 2 = y 1 ? x 1 x_2=y_1-x_1 x2?=y1??x1?，可以省略它)：
E X h [ L L ( θ ∣ X , Y ) ] = E X h [ ∑ i = 1 N l nP ( x i 1 , x i 2 , x i 3 ∣ y i 1 , y i 2 ; θ ) ] = ∑ i = 1 N E X h [ l nP ( x i 1 , x i 2 , x i 3 ∣ y i 1 , y i 2 ; θ ) ] \begin{aligned}\mathbb{E}_{X_h}[LL(\pmb\theta|X,Y)] & = \mathbb{E}_{X_h}[\sum\limits_{i=1}^{N}ln\ P(x_{i1},x_{i2},x_{i3}|y_{i1},y_{i2}; \theta)] \\ & = \sum\limits_{i=1}^{N}\mathbb{E}_{X_h}[ln\ P(x_{i1},x_{i2},x_{i3}|y_{i1},y_{i2}; \theta)]\end{aligned} EXh??[LL(θθ∣X,Y)]?=EXh??[i=1∑N?ln P(xi1?,xi2?,xi3?∣yi1?,yi2?; θ)]=i=1∑N?EXh??[ln P(xi1?,xi2?,xi3?∣yi1?,yi2?; θ)]?
其中
E X h [ l nP ( x i 1 , x i 2 , x i 3 ∣ y i 1 , y i 2 ; θ ) ] = E X h [ l n( x 1 + x 2 + x 3 ) ! x 1 ! x 2 ! x 3 ! ( 1 4 ) x 1 ( 1 + θ 4 ) x 2 ( 2 ? θ 4 ) x 3 ] = E X h [ l n( x 1 + x 2 + x 3 ) ! x 1 ! x 2 ! x 3 ! + x 1 l n 1 4 + x 2 l n ( 1 + θ 4 ) + x 3 l n ( 2 ? θ 4 ) ] = E X h [ l n( x 1 + x 2 + x 3 ) ! x 1 ! x 2 ! x 3 ! ] + l n 1 4 E X h [ x 1 ] + l n ( 1 + θ 4 ) E X h [ x 2 ] + l n ( 2 ? θ 4 ) E X h [ x 3 ] . \begin{aligned}\mathbb{E}_{X_h}[ln\ P(x_{i1},x_{i2},x_{i3}|y_{i1},y_{i2}; \theta)] & = \mathbb{E}_{X_h}[ln \ \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}(\frac14)^{x_1}(\frac{1+\theta}{4})^{x_2}(\frac{2-\theta}{4})^{x_3}]\\ & =\mathbb{E}_{X_h}[ln \ \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}+x_1ln\frac14+x_2ln(\frac{1+\theta}{4})+x_3ln(\frac{2-\theta}{4})]\\ & = \mathbb{E}_{X_h}[ln \ \frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}] +ln\frac14\mathbb{E}_{X_h}[x_1]+ln(\frac{1+\theta}{4})\mathbb{E}_{X_h}[x_2]+ln(\frac{2-\theta}{4})\mathbb{E}_{X_h}[x_3]\end{aligned}. EXh??[ln P(xi1?,xi2?,xi3?∣yi1?,yi2?; θ)]?=EXh??[ln x1?!x2?!x3?!(x1?+x2?+x3?)!?(41?)x1?(41+θ?)x2?(42?θ?)x3?]=EXh??[ln x1?!x2?!x3?!(x1?+x2?+x3?)!?+x1?ln41?+x2?ln(41+θ?)+x3?ln(42?θ?)]=EXh??[ln x1?!x2?!x3?!(x1?+x2?+x3?)!?]+ln41?EXh??[x1?]+ln(41+θ?)EXh??[x2?]+ln(42?θ?)EXh??[x3?]?.
（3）求得使(2)中期望最大的 θ \theta θ值,为此需要对期望求 θ \theta θ的导数：
∑ i = 1 N d E X h d θ = ∑ i = 1 N ( 1 θ + 1 E X h [ x i 2 ] ? 1 2 ? θ E X h [ x i 3 ] ) = 0 \begin{aligned}\sum\limits_{i=1}^{N}\frac{d\mathbb{E}_{X_h}}{d\theta} & = \sum\limits_{i=1}^{N}(\frac{1}{\theta+1}\mathbb{E}_{X_h}[x_{i2}]-\frac{1}{2-\theta}\mathbb{E}_{X_h}[x_{i3}])=0\\ \end{aligned} i=1∑N?dθdEXh????=i=1∑N?(θ+11?EXh??[xi2?]?2?θ1?EXh??[xi3?])=0?
得
θ t + 1 = ∑ i = 1 N ( 2 E X h [ x i 2 ] ? E X h [ x i 3 ] ) ∑ i = 1 N ( E X h [ x i 2 ] + E X h [ x i 3 ] ) . \theta^{t+1}=\frac{\sum\limits_{i=1}^{N}(2\mathbb{E}_{X_h}[x_{i2}]-\mathbb{E}_{X_h}[x_{i3}])}{\sum\limits_{i=1}^{N}(\mathbb{E}_{X_h}[x_{i2}]+\mathbb{E}_{X_h}[x_{i3}])}. θt+1=i=1∑N?(EXh??[xi2?]+EXh??[xi3?])i=1∑N?(2EXh??[xi2?]?EXh??[xi3?])?.
由于 E X h [ ? ] \mathbb{E}_{X_h}[\cdot] EXh??[?]表示已知 θ t , y i 1 , y i 2 \theta^t,y_{i1},y_{i2} θt,yi1?,yi2?的对 x 1 x_1 x1?的期望，因此
E X h [ x i 3 ∣ y i 1 , y i 2 ] = y i 2 ( x i 3 = y i 2 ) ; E X h [ x i 2 ∣ y i 1 , y i 2 ] = y i 1 1 + θ t 4 1 4 + 1 + θ t 4 = y i 1 1 + θ t 2 + θ t ( 见下方说明 ) . \mathbb{E}_{X_h}[x_{i3}|y_{i1},y_{i2}]=y_{i2}(x_{i3}=y_{i2}); \\ \\ \mathbb{E}_{X_h}[x_{i2}|y_{i1},y_{i2}]=y_{i1}\frac{\frac{1+\theta^t}{4}}{\frac14+\frac{1+\theta^t}{4}}=y_{i1}\frac{1+\theta^t}{2+\theta^t}(见下方说明). EXh??[xi3?∣yi1?,yi2?]=yi2?(xi3?=yi2?); EXh??[xi2?∣yi1?,yi2?]=yi1?41?+41+θt?41+θt??=yi1?2+θt1+θt?(见下方说明).
对于上式的第二行，原论文的第53页附加了一个简单的证明。这里我们不作严格的证明，仅从三项分布的意义和性质来理解：三项分布的一般形式如
P ( x 1 , x 2 , x 3 ) = ( x 1 + x 2 + x 3 ) ! x 1 ! x 2 ! x 3 ! p x 1 q x 2 r x 3 P(x_1,x_2,x_3)=\frac{(x_1+x_2+x_3)!}{x_1!x_2!x_3!}p^{x_1}q^{x_2}r^{x_3} P(x1?,x2?,x3?)=x1?!x2?!x3?!(x1?+x2?+x3?)!?px1?qx2?rx3?
其中 p + q + r = 1 , p+q+r=1, p+q+r=1,分别代表在总体中抽出第一/二/三类样本的概率为 p / q / r p/q/r p/q/r（放回抽样，和二项分布很类似）；整个式子代表抽取 n = x 1 + x 2 + x 3 n=x_1+x_2+x_3 n=x1?+x2?+x3?次，分别抽得一/二/三类样本 x 1 / x 2 / x 3 x_1/x_2/x_3 x1?/x2?/x3?个的概率。回到上面第二行的问题：当 x i 1 + x i 2 x_{i1}+x_{i2} xi1?+xi2?固定为 y i 1 y_{i1} yi1?时,我们想求得 E X h [ x i 2 ∣ y i 1 , y i 2 ] . \mathbb{E}_{X_h}[x_{i2}|y_{i1},y_{i2}]. EXh??[xi2?∣yi1?,yi2?].这时相当于一个总数为 y i 1 , y_{i1}, yi1?,概率为 1 + θ t 4 1 4 + 1 + θ t 4 \frac{\frac{1+\theta^t}{4}}{\frac14+\frac{1+\theta^t}{4}} 41?+41+θt?41+θt??(将 x 1 x_1 x1?和 x 2 x_2 x2?的概率进行“加权”分配，即归一化)的二项分布，即给定 y i 1 y_{i1} yi1?的条件下 x i 2 ～ B ( y i 1 , 1 + θ t 2 + θ t ) . x_{i2}\sim B(y_{i1},\frac{1+\theta^t}{2+\theta^t}). xi2?～B(yi1?,2+θt1+θt?).由二项分布的性质得
E X h [ x i 2 ∣ y i 1 , y i 2 ] = y i 1 1 + θ t 2 + θ t . \mathbb{E}_{X_h}[x_{i2}|y_{i1},y_{i2}]=y_{i1}\frac{1+\theta^t}{2+\theta^t}. EXh??[xi2?∣yi1?,yi2?]=yi1?2+θt1+θt?.
综上所述，参数更新表达式为
θ t + 1 = ∑ i = 1 N ( 2 y i 1 1 + θ t 2 + θ t ? y i 2 ) ∑ i = 1 N ( y i 1 1 + θ t 2 + θ t + y i 2 ) . \theta^{t+1}=\frac{\sum\limits_{i=1}^{N}(2y_{i1}\frac{1+\theta^t}{2+\theta^t}-y_{i2})}{\sum\limits_{i=1}^{N}(y_{i1}\frac{1+\theta^t}{2+\theta^t}+y_{i2})}. θt+1=i=1∑N?(yi1?2+θt1+θt?+yi2?)i=1∑N?(2yi1?2+θt1+θt??yi2?)?.
下面我们也是用python来验证一下这个结果：

import random import numpy as np# 每次三项分布总共取出20个样本 n = 20 # 参数theta theta = 1.5# 返回观测值y_1, y_2 def observe(n, p1, p2, p3): sum1, sum2, sum3 = 0, 0, 0 for i in range(n): t = random.random() if t < p1: sum1 += 1 elif t < p1 + p2: sum2 += 1 else: sum3 += 1 return [sum1 + sum2, sum3]# N 为数据集大小 def get_dataset(N): dataset = [] for i in range(N): dataset.append(observe(n, 1 / 4, (1 + theta) / 4, (2 - theta) / 4)) return datasetdef EM(dataset): # theta在(-1, 2)之间 theta_hat = random.random() * 3 - 1 y = np.array(dataset, dtype = 'float')# A,B 分别为y_i1,y_i2之和 A = np.sum(y[:, 0]) B = np.sum(y[:, 1])iteration = 0 while True: iteration += 1tmp = (1 + theta_hat) / (2 + theta_hat) new_theta_hat = (2 * A * tmp - B) / (A * tmp + B)conv = abs(new_theta_hat - theta_hat) if conv < 1e-10 or iteration > 1000: breaktheta_hat = new_theta_hatprint(f'iteration = {iteration}') return theta_hatdataset = get_dataset(10000) for i in range(20): print(EM(dataset))

运行结果：

iteration = 8 1.5024600000467414 iteration = 9 1.5024599999731247 iteration = 9 1.5024600000072472 iteration = 10 1.502459999983291 iteration = 9 1.5024600000063015 # ...(后面省略)

结果普遍收敛到初始定下的 θ = 1.5 , \theta=1.5, θ=1.5,且收敛速度很快。
参考文献

Do, C., Batzoglou, S. What is the expectation maximization algorithm?. Nat Biotechnol 26, 897–899 (2008). ??
A. O. Hero, “On the Convergence of the Em Algorithm,” Proceedings. IEEE International Symposium on Information Theory, 1993, pp. 187-187, doi: 10.1109/ISIT.1993.748501. ??
【算法|计算建模之EM算法】Molenberghs G , Verbeke G . The Expectation-Maximization Algorithm[M]. 2000. ??