《机器学习公式推导与代码实现》chapter22-EM算法

《机器学习公式推导与代码实现》学习笔记，记录一下自己的学习过程，详细的内容请大家购买作者的书籍查阅。

EM算法

作为一种迭代算法，EM算法(expectation maximization，期望极大值算法)用于包含隐变量的概率模型参数的极大似然估计。

EM算法包括两个步骤：E步，求期望(expectation)；M步，求极大(maximization)。

1 极大似然估计

极大似然估计（maximum likelihood estimation,MLE）是统计学领域中一种经典的参数估计方法。对于某个随机样本满足某种概率分布，但其中的统计参数未知的情况，极大似然估计可以让我通过若干次试验的结果来估计参数的值。

以一个经典的例子进行说明，比如我们想了解某高校学生的身高分布。我们假设该校学生的身高分布服从一个正态分布$N(\mu ,{\sigma }^2)$,其中分布参数$\mu$和${\sigma }^2$未知。全校有数万名学生，要一个一个实测肯定不现实，所以我们决定采用统计抽样的方法，随机选取100名学生测得其身高。

要通过100人的身高估算全校学生的身高，需要明确以下问题。第一个问题是抽到这100人的概率是多少。因为每个人的选取都是独立的，所以抽到这100人的概率可以表示为单个概率的乘积：
$$L(\theta )=L(x_1,x_2,\cdots ,x_n;\theta )=\prod _{i=1}^{n}p(x_i\mid \theta )$$
上式为似然函数，为了计算方便，我们会对似然函数取对数：
$$H(\theta )={\ln }_{}L(\theta )={\ln }_{}\prod _{i=1}^{n}p(x_i\mid \theta )=\sum _{i=1}^n{\ln }_{}p(x_i\mid \theta )$$
第二个问题是为什么刚好抽到这100人。按照极大似然估计理论，在学校这么多学生中，我们恰好抽到这100人而不是另外100人，正是因为这100人出现的概率极大，即其对应的似然函数极大：
$$\hat{\theta }=argmaxL(\theta )$$
最后一个问题是如何求解，直接对$L(\theta )$求导并使其为0。

所以极大似然估计法可以看作由抽样结果对条件的反推，即已知某个参数能使得这些样本出现的概率极大，我们就直接把该参数作为参数估计的真实值。

2 EM算法

假设全校学生的身高付出一个分布的假设过于笼统，实际上男女分布不同，假设其中男生身高服从分布$N({\mu }_1^{},{\sigma }_1^2)$,女生身高分布为$N({\mu }_2^{},{\sigma }_2^2)$。现在估计该校学生身高，就不能简单地使用一个分布的假设了。

假设分别抽取50个男生和50个女生，对他们分开进行估计。假设我们并不知道抽样得到的这样样本来自男生还是女生。

学生的身高是观测变量(observable variable)，样本的性别是一种隐变量(hidden variable)。

现在我们需要估计两个问题：一是这个样本是男生的还是女生的，而是男生和女生对应的身高的正态分布参数分别是多少。这种情况极大似然估计就不太适用了，要估计男女生身高分布，就必须先估计该学生是男是女。反过来要估计该学生是男还是女，又得从身高来判断。但二者相互依赖，直接用极大似然估计无法计算。

针对这种包含隐变量的参数估计问题，一般使用EM(expectation maximization)算法，即期望极大化算法来进行求解。针对上述身高估计问题，EM算法的求解思路是：既然两个问题相互依赖，这肯定是一个动态求解过程。不如直接给定男女身高的分布初始值，根据初始值估计哪个样本是男/女生的概率(E步)，然后据此使用极大似然估计男女生的身高分布参数(M步)，之后动态迭代调整到满足终止条件为止。

EM算法的应用场景就是解决包含隐变量的概率模型参数估计问题。给定观测变量数据Y，隐变量数据Z，联合概率分布$P(Y,Z|\theta )$，以及关于隐变量的条件分布$P(Z|Y,\theta )$，使用EM算法对模型参数$\theta$进行估计的流程如下：
(1)初始化模型参数${\theta }^{(0)}$，开始迭代。
(2)E步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算Q函数：
$$Q(\theta ,{\theta }^{(i)})=E_Z\left[{\log }_{}P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]=\sum _Z{\log }_{}P(Y,Z\mid \theta )P(Z|Y,{\theta }^{(i)})$$
其中$P(Z|Y,{\theta }^{(i)})$为给定观测数据$Y$和当前参数估计${\theta }^{(i)}$的情况下隐变量数据$Z$的条件概率分布。E步的关键是这个Q函数，Q函数定义为完全数据的对数似然函数${\log }_{}P(Y,Z\mid \theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$的情况下未观测数据$Z$的条件概率分布。
(3)M步：求使得Q函数最大化的参数$\theta$，确定第$i+1$次迭代的参数估计值${\theta }^{(i+1)}$：
$${\theta }^{(i+1)}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{(i)})$$
(4)重复迭代E步和M步直至收敛。

由EM算法过程可知，其关键在于E步要确定Q函数。E步在固定模型参数的情况下估计隐状态变量分布，而M步则是固定隐变量来估计模型参数。二者交互进行，直至算法收敛条件。

EM算法动态迭代过程：

3 三硬币模型

4 基于numpy实现三硬币模型

import numpy as np

## EM算法过程定义
def em(data, thetas, max_iter=50, eps=1e-3): # data观测数据，thetas初始化的估计参数值，eps收敛阈值

    ll_old = 0 # 初始化似然函数值
    for i in range(max_iter):
        # E步：求隐变量分布
        log_like = np.array([np.sum(data*np.log(theta), axis=1) for theta in thetas]) # 对数似然 2*5
        like = np.exp(log_like) # 似然 2*5
        ws = like/like.sum(0) # 隐变量分布 2*5
        ll_new = np.sum([w*l for w, l in zip(ws, log_like)]) # 期望

        # M步：更新参数值
        vs = np.array([w[:, None] * data for w in ws]) # 概率加权 2*5*2
        thetas = np.array([v.sum(0)/v.sum() for v in vs]) # 2*2

        # 打印结果
        print(f'Iteration:{i+1}')
        print(f'theta_B = {thetas[0,0]:.2}, theta_C = {thetas[1,0]:.2}, {ll_new}')

        # 满足条件退出迭代
        if np.abs(ll_new - ll_old) < eps:
            break
        ll_old = ll_new
    
    return thetas

EM算法求解三硬币问题：

# 观测数据，5次独立实验，每次试验10次抛掷的正反面次数
observed_data = np.array([(5, 5), (9, 1), (8, 2), (4, 6), (7, 3)]) # 比如第一次试验为5次正面5次反面
# 初始化参数值，硬币B出现正面的概率0.6，硬币C出现正面的概率为0.5
thetas = np.array([[0.6, 0.4], [0.5, 0.5]])
# EM算法寻优
thetas = em(observed_data, thetas, max_iter=30)
thetas

Iteration:1
theta_B = 0.71, theta_C = 0.58, -32.68721052517165
Iteration:2
theta_B = 0.75, theta_C = 0.57, -31.258877917413145
Iteration:3
theta_B = 0.77, theta_C = 0.55, -30.760072598843628
Iteration:4
theta_B = 0.78, theta_C = 0.53, -30.33053606687176
Iteration:5
theta_B = 0.79, theta_C = 0.53, -30.071062062760774
Iteration:6
theta_B = 0.79, theta_C = 0.52, -29.95042921516964
Iteration:7
theta_B = 0.8, theta_C = 0.52, -29.90079955867412
Iteration:8
theta_B = 0.8, theta_C = 0.52, -29.881202814860167
Iteration:9
theta_B = 0.8, theta_C = 0.52, -29.873553692091832
Iteration:10
theta_B = 0.8, theta_C = 0.52, -29.870576075992844
Iteration:11
theta_B = 0.8, theta_C = 0.52, -29.86941691676721
Iteration:12
theta_B = 0.8, theta_C = 0.52, -29.868965223428773

array([[0.7967829 , 0.2032171 ],
       [0.51959543, 0.48040457]])

算法在第7次迭代时收敛，最后硬币B和硬币C出现正面的概率分别为0.80和0.52。

笔记本_Github地址