09. EM

算法思路

之前所遇到的模型，大多是求$\hat{\theta }=\arg {\max }_{\theta }P(X,y|\theta )$，但是在一些存在隐变量的模型中，表达式$P(X,y|\theta )$难以求解析解，所以可以通过一个迭代的算法求出其解，这个迭代的算法就称为EM算法。

算法流程

1. 确定初始参数${\theta }^{(0)}$
2. 计算$Q(\theta ,{\theta }^{(i)})$
   $$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& ={\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}[\log P(X,Z|\theta )]\end{array}$$
3. 求解${\theta }^{(i+1)}=\arg {\max }_{\theta }Q(\theta ,{\theta }^{(i)})$

收敛性证明

我们只需要证明$P(X|{\theta }^{(i+1)})>P(X|{\theta }^{(i)})$即可。
$$\begin{array}{rl}P(Z|X,\theta )\ast P(X|\theta )& =P(X,Z|\theta )\\ \log P(X|\theta )& =\log P(X,Z|\theta )-\log P(Z|X,\theta )\\ {\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(X|\theta )& ={\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(X,Z|\theta )-{\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(Z|X,\theta )\\ \log P(X|\theta )& ={\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(X,Z|\theta )-{\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(Z|X,\theta )\\ \log P(X|\theta )& =Q(\theta ,{\theta }^{(i)})-{\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log P(Z|X,\theta )\end{array}$$
所以$D\_value=\log P(X|{\theta }^{(i+1)})-\log P(X|{\theta }^{(i)})$的值为：
$$\begin{array}{rl}D\_value& =Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})+{\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log \frac{P(Z|X,{\theta }^{(i)})}{P(Z|X,{\theta }^{(i+1)})}\\ & =Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})+D_{KL}(P(Z|X,{\theta }^{(i)})\Vert P(Z|X,{\theta }^{(i+1)}))\end{array}$$
由于$Q({\theta }^{(i+1)},{\theta }^{(i)})=\max Q(\theta ,{\theta }^{(i)})\ge Q({\theta }^{(i)},{\theta }^{(i)})$，又因为差值$D\_value\ge 0$，所以算法收敛。

算法导出

虽然我们知道EM算法流程，但是我们不知道如何导出EM算法。
仿照贝叶斯神经网络时的推导，我们假设$Z\sim q(Z)$，那么我们对等式$\log P(X|\theta )=\log P(X,Z|\theta )-\log P(Z|X,\theta )$求期望得：
$$\begin{array}{rl}\log P(X|\theta )& =\log P(X,Z|\theta )-\log P(Z|X,\theta )\\ & =\log \frac{P(X,Z|\theta )}{q(Z)}-\log \frac{P(Z|X,\theta )}{q(Z)}\\ & ={\mathbb{E}}_{q(Z)}\log \frac{P(X,Z|\theta )}{q(Z)}-D_{KL}(q(Z)\Vert P(Z|X,\theta ))\end{array}$$
我们知道当$q(Z)=P(Z|X,\theta )$时，$D_{KL}=0$。令$q(Z)=P(Z|X,\theta )$得：
$$\begin{array}{rl}\log P(X|\theta )& ={\mathbb{E}}_{P(Z|X,\theta )}\log \frac{P(X,Z|\theta )}{P(Z|X,\theta )}\end{array}$$
但是这一步仍然不可解，所以采用迭代算法策略：$\log P(X|\theta )={\mathbb{E}}_{P(Z|X,{\theta }^{(i)})}\log \frac{P(X,Z|\theta )}{P(Z|X,{\theta }^{(i)})}$，由于$\log$下面得分母与优化过程无关，所以略去。
所以导出EM算法。

代码实现

import numpy as np
import random

def load_data(mus, sigmas, alphas, length):
    assert(len(mus) == len(sigmas) and len(mus) == len(alphas))
    lengths = [round(length * alpha) for alpha in alphas]
    data = []
    for mu, sigma, length in zip(mus, sigmas, lengths):
        data.extend(np.random.normal(mu, sigma, length))
    random.shuffle(data)
    return np.array(data)


def gaussian(data, mu, sigma):
    return np.exp(-(data - mu)**2 / 2 / sigma**2) / np.sqrt(2 * np.pi) / sigma


def train(data, mus, sigmas, alphas, epochs):
    for epoch in range(epochs):
        gammas, gamma_sum = [], 0
        for alpha, mu, sigma in zip(alphas, mus, sigmas):
            gamma = alpha * gaussian(data, mu, sigma)
            gammas.append(gamma)
            gamma_sum += gamma
        for i in range(len(mus)):
            gammas[i] /= gamma_sum
            mus[i] = np.dot(gammas[i], data) / np.sum(gammas[i])
            sigmas[i] = np.sqrt(np.dot(gammas[i], (data - mus[i])**2) / np.sum(gammas[i]))
            alphas[i] = np.sum(gammas[i]) / len(gammas[i])
    return mus, sigmas, alphas


if __name__ == '__main__':
    mus, sigmas, alphas = [0.7, 0.5], [1, 2], [0.3, 0.7]
    data = load_data(mus, sigmas, alphas, 10000)
    mus, sigmas, alphas = train(data, [0, 1], [1, 1], [0.5, 0.5], 10000)
    print(mus, sigmas, alphas)