期望极大（EM）算法

EM算法

概率模型有时既含有观测变量（observable vriable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。

EM算法[Dempster, 1977]是一种迭代算法，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximization）。所以这一算法称为期望极大算法（expectation maximization algorithm，EM）

1 EM 算法

import numpy as np

def e_step(y, pi, p, q):
    
    mu_1 = pi * p ** y * (1 - p) ** (1 - y)
    mu_2 = (1 - pi) * q ** y * (1 - q) ** (1 - y)
    
    mu = mu_1 / (mu_1 + mu_2)
    
    return mu

def m_step(y, mu):
    
    n = len(y)
    pi = np.sum(mu) / n
    p = sum(y * mu) / sum(mu)
    q = sum(y * (1 - mu)) / sum(1 - mu)
    
    return pi, p, q

def diff(pi, p, q, pi_, p_, q_):
    
    return np.sum(np.abs([pi - pi_, p - p_, q - q_]))

def em(y, pi, p, q):
    cnt = 1
    while True:

        print("-" * 10)
        print("iter %d:" % cnt)
        pi_ = pi
        p_ = p
        q_ = q

        mu = e_step(y, pi, p, q)
        print(mu)
        pi, p, q = m_step(y, mu)
        print(pi, p, q)

        if diff(pi, p, q, pi_, p_, q_) < 0.001:
            break

        cnt += 1
        
    return pi, p, q

y = np.array([1, 1, 0, 1, 0, 0, 1, 0, 1, 1])

print("*" * 10)
pi = 0.5
p = 0.5
q = 0.5

pi, p, q = em(y, pi, p, q)

print("*" * 10)
pi = 0.4
p = 0.6
q = 0.7

pi, p, q = em(y, pi, p, q)

print("*" * 10)
pi = 0.46
p = 0.55
q = 0.67

pi, p, q = em(y, pi, p, q)

**********
----------
iter 1:
[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
0.5 0.6 0.6
----------
iter 2:
[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]
0.5 0.6 0.6
**********
----------
iter 1:
[0.36363636 0.36363636 0.47058824 0.36363636 0.47058824 0.47058824
 0.36363636 0.47058824 0.36363636 0.36363636]
0.40641711229946526 0.5368421052631579 0.6432432432432431
----------
iter 2:
[0.36363636 0.36363636 0.47058824 0.36363636 0.47058824 0.47058824
 0.36363636 0.47058824 0.36363636 0.36363636]
0.40641711229946526 0.5368421052631579 0.6432432432432431
**********
----------
iter 1:
[0.41151594 0.41151594 0.53738318 0.41151594 0.53738318 0.53738318
 0.41151594 0.53738318 0.41151594 0.41151594]
0.461862835113919 0.5345950037850112 0.6561346417857326
----------
iter 2:
[0.41151594 0.41151594 0.53738318 0.41151594 0.53738318 0.53738318
 0.41151594 0.53738318 0.41151594 0.41151594]
0.46186283511391907 0.5345950037850112 0.6561346417857326

通常情况，$Y$表示观测随机变量的数据，$Z$表示隐随机变量的数据。$Y$和$Z$均已知称为完全数据（complete-data），仅有观测数据$Y$称为不完全数据（incomplete-data）。假设给定观测数据$Y$，其概率分布是$P(Y;\theta )$，其中$\theta$是需要估计的模型参数，那么不完全数据$Y$的似然函数是$P(Y;\theta )$，对数似然函数$L(\theta )=\log P(Y;\theta )$；假设$Y$和$Z$的联合概率分布是$P(Y,Z;\theta )$，那么完全数据的对数似然函数是$L(\theta )=\log P(Y,Z;\theta )$。

EM算法通过迭代求解$L(\theta )=\log P(Y,Z;\theta )$的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。

算法9.1（EM算法）

输入：观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z;\theta )$，条件分布$P(Z|Y;\theta )$；

输出：模型参数$\theta$。

1. 选择参数初值${\theta }^{(0)}$，开始迭代；

2. E步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算

$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& ={\text{E}}_Z\left[\log P(Y,Z;\theta )|Y;{\theta }^{(i)}\right]\\ & =\sum _{Z}P(Z|Y;{\theta }^{(i)})\log P(Y,Z;\theta )\end{array}\text{\hspace{1em}(9)}$$

其中，$P(Z;Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布；

3. M步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$

$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(10)}$$

4. 重复第2步和第3步，直到收敛。

方程（9）的函数$Q(\theta ,{\theta }^{(i)})$是EM算法的核心，称为$Q$函数（$Q$ function）。

定义9.1（Q函数） 完全数据的对数似然函数$P(Y,Z|\theta )$是关于给定观测数据$Y$和当前参数${\theta }^{(i)}$下，对未观测数据$Z$条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望，称为$Q$函数，即

Q ( θ , θ ( i ) ) = E Z [ log ⁡ P ( Y , Z ; θ ) ∣ Y ; θ