EM算法及python简单实现

最大期望算法（Expectation-maximization algorithm，又译为期望最大化算法），是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐性变量z。

算法为：
0.初始化概率分布参数
1.估计隐变量分布

2.最大化下式

3.循环上述两步

举例：给定一组数据。每个数据可能来自于两种正态分布，其中每种正态分布均值不知道，方差已知。
以下是代码：

import numpy as np

# 数据
x1 = np.random.normal(loc=20, size=[1000])
x2 = np.random.normal(loc=-10, size=[100])
x = np.hstack([x1, x2])

# 保存后验概率
p1 = np.zeros_like(x)
p2 = p1.copy()

# 均值初始值
mu1 = -1
mu2 = 1
#先验概率初始值
pi1, pi2 = 0.1, 0.9


def getp(x, mu, sigma):
    return 1 / (2 * np.pi) ** (1 / 2) / sigma * np.exp(-1 / 2 * (x - mu) ** 2 / 2 / sigma ** 2)


for i in range(100):
    p1 = getp(x, mu1, sigma=1)
    p2 = getp(x, mu2, sigma=2)

    # E步，求得隐变量后验概率
    p = p1 / (p1 + p2)
    pi1 = sum(p) / len(p)
    pi2 = 1 - pi1

    # M步，极大似然估计参数
    mu1 = (p * x).sum() / sum(p)
    mu2 = ((1 - p) * x).sum() / sum(1 - p)

print(mu1, mu2)
print(pi1, pi2)