EM算法

EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法

含有隐变量的概率模型的数据表示为θ\thetaθ

这里，YYY是观测变量的数据，ZZZ是隐变量的数据，θ\thetaθ 是模型参数

EM算法通过迭代求解观测数据的对数似然函数

L(θ)=log⁡P(Y∣θ){L}(\theta)=\log {P}(\mathrm{Y} | \theta)L(θ)=logP(Y∣θ)的极大化，实现极大似然估计

每次迭代包括两步

EEE步，求期望

即求logP(Z∣Y,θ)logP\left(Z | Y, \theta\right)logP(Z∣Y,θ) )关于P(Z∣Y,θ(i))P\left(Z | Y, \theta^{(i)}\right)P(Z∣Y,θ(i)))的期望：

Q(θ,θ(i))=∑Zlog⁡P(Y,Z∣θ)P(Z∣Y,θ(i))Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right)Q(θ,θ(i))=Z∑​logP(Y,Z∣θ)P(Z∣Y,θ(i))

称为QQQ函数

这里θ(i)\theta^{(i)}θ(i)是参数的现估计值

MMM步，求极大

即极大化QQQ函数得到参数的新估计值

θ(i+1)=arg⁡max⁡θQ(θ,θ(i))\theta^{(i+1)}=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right)θ(i+1)=argθmax​Q(θ,θ(i))

在构建具体的EM算法时，重要的是定义QQQ函数

每次迭代中，EM算法通过极大化QQQ函数来增大对数似然函数L(θ){L}(\theta)L(θ)

EMEMEM算法在每次迭代后均提高观测数据的似然函数值，即

P(Y∣θ(i+1))⩾P(Y∣θ(i))P\left(Y | \theta^{(i+1)}\right) \geqslant P\left(Y | \theta^{(i)}\right)P(Y∣θ(i+1))⩾P(Y∣θ(i))

在一般条件下EM算法是收敛的，但不能保证收敛到全局最优。

EMEMEM算法应用极其广泛

主要应用于含有隐变量的概率模型的学习

高斯混合模型的参数估计是EM算法的一个重要应用

隐马尔可夫模型的非监督学习也是EM算法的一个重要应用

EMEMEM算法还可以解释为FFF函数的极大-极大算法

EMEMEM算法有许多变形，如GEM算法

GEMGEMGEM算法的特点是每次迭代增加FFF函数值（并不一定是极大化FFF函数），从而增加似然函数值

P(Y∣θ)=∏[πpyi(1−p)1−yi+(1−π)qyi(1−q)1−yi]P(Y|\theta) = \prod[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi) q^{y_i}(1-q)^{1-y_i}]P(Y∣θ)=∏[πpyi​(1−p)1−yi​+(1−π)qyi​(1−q)1−yi​]

代码实现

E step:

μi+1=π(pi)yi(1−(pi))1−yiπ(pi)yi(1−(pi))1−yi+(1−π)(qi)yi(1−(qi))1−yi\mu^{i+1}=\frac{\pi (p^i)^{y_i}(1-(p^i))^{1-y_i}}{\pi (p^i)^{y_i}(1-(p^i))^{1-y_i}+(1-\pi) (q^i)^{y_i}(1-(q^i))^{1-y_i}}μi+1=π(pi)yi​(1−(pi))1−yi​+(1−π)(qi)yi​(1−(qi))1−yi​π(pi)yi​(1−(pi))1−yi​​

def pmf(i, pro_A, pro_B, por_C):pro_1 = pro_A * math.pow(pro_B, data[i]) * math.pow((1 - pro_B), 1 - data[i])pro_2 = pro_A * math.pow(pro_C, data[i]) * math.pow((1 - pro_C), 1 - data[i])return pro_1 / (pro_1 + pro_2)

M step:

πi+1=1n∑j=1nμji+1\pi^{i+1}=\frac{1}{n}\sum_{j=1}^n\mu^{i+1}_jπi+1=n1​j=1∑n​μji+1​

pi+1=∑j=1nμji+1yi∑j=1nμji+1p^{i+1}=\frac{\sum_{j=1}^n\mu^{i+1}_jy_i}{\sum_{j=1}^n\mu^{i+1}_j}pi+1=∑j=1n​μji+1​∑j=1n​μji+1​yi​​

qi+1=∑j=1n(1−μji+1yi)∑j=1n(1−μji+1)q^{i+1}=\frac{\sum_{j=1}^n(1-\mu^{i+1}_jy_i)}{\sum_{j=1}^n(1-\mu^{i+1}_j)}qi+1=∑j=1n​(1−μji+1​)∑j=1n​(1−μji+1​yi​)​

class EM:def __init__(self, prob):self.pro_A, self.pro_B, self.pro_C = prob# e_stepdef pmf(self, i):pro_1 = self.pro_A * math.pow(self.pro_B, data[i]) * math.pow((1 - self.pro_B), 1 - data[i])pro_2 = (1 - self.pro_A) * math.pow(self.pro_C, data[i]) * math.pow((1 - self.pro_C), 1 - data[i])return pro_1 / (pro_1 + pro_2)# m_stepdef fit(self, data):count = len(data)print('init prob:{}, {}, {}'.format(self.pro_A, self.pro_B,self.pro_C))for d in range(count):_ = yield_pmf = [self.pmf(k) for k in range(count)]pro_A = 1 / count * sum(_pmf)pro_B = sum([_pmf[k] * data[k] for k in range(count)]) / sum([_pmf[k] for k in range(count)])pro_C = sum([(1 - _pmf[k]) * data[k]for k in range(count)]) / sum([(1 - _pmf[k])for k in range(count)])print('{}/{}  pro_a:{:.3f}, pro_b:{:.3f}, pro_c:{:.3f}'.format(d + 1, count, pro_A, pro_B, pro_C))self.pro_A = pro_Aself.pro_B = pro_Bself.pro_C = pro_C