EM算法解析与实践-优快云博客

本文链接：https://blog.youkuaiyun.com/yy2050645/article/details/82596704

EM算法是一种迭代求解隐变量概率模型参数的一种算法，其本质就是使用极大似然估计，但与以往不同的是，它并没有直接对似然函数求极大，而是以求取似然函数的下界的极大来近似逼近似然函数的极大，那么接下来就详细介绍一下这个算法。

在介绍这个算法之前，首先介绍一下Jensen不等式：

首先我们给定一个函数 $f(x)$ 形状如下所示：

显然它是一个凸函数，那么凸函数的定义是什么呢？如图所示，我们可以得凸函数的定义就是：

$f(\alpha*a+(1-\alpha) * b)\leq \alpha * f(a)+(1-\alpha)*f(b)$
由于 $0<\alpha<1$ ，因此我们可以将 $\alpha$ 和 $1-\alpha$ 看做是一个概率分布，在不等号左边，可以将其看作是 $f(x)$ 中，随机变量 $x$ 的分布，此时随机变量只能取两个值 $a$ 或 $b$ ，如下所示：

$x$	$a$	$b$
$p(x)$	$\alpha$	$1-\alpha$

同理，在不等号右边，可以将其看作是随机变量 $f(x)$ 的概率分布， $f(x)$ 只能取值 $f(a)$ 以及 $f(b)$ ，如下所示：

$f(x)$	$f(a)$	$f(b)$
$p(f(x))$	$\alpha$	$1-\alpha$

那么，凸函数定义表达式便可写为：

$f(E(x))\leq E(f(x))$

式中 $E$ 表示期望，这便是Jensen不等式，也就是说，如果一个函数 $f(x)$ 为凸函数，则其满足上式。

同理我们可以推出如果一个函数 $f(x)$ 为凹函数，则满足下式：

$f(E(x))\geq E(f(x))$

Jensen不等式便介绍完毕了，我们接着介绍EM算法。

有时概率模型既含有观测变量，也含有隐变量，隐变量就是不可观测的变量，其值我们是观测不到的，那么如何来求解隐变量的概率模型的参数呢？

下面我们先来给出一个含有隐变量的场景：

三硬币模型：假设有三枚硬币，分别记做 $A，B，C$ ， $B$ ， $C$ 。当我们抛掷这些硬币时，其正面朝上的概率依次为 $\pi$ ， $p$ 以及 $q$ 。接下来我们进行如下实验：第一轮我们先抛掷硬币 $A，B，C$ ，之后根据硬币 $A，B，C$ 的结果我们选择第二轮抛掷硬币 $B$ 还是硬币 $C$ 。如果第一轮硬币 $A，B，C$ 出现的是正面，我们第二轮选择抛掷硬币 $B$ ；如果第一轮硬币 $A，B，C$ 出现的是反面，第二轮我们就抛掷硬币 $C$ 。如果抛掷一枚硬币出现正面朝上我们记做1，反面朝上记做0。假设我么独立的重复进行10次试验，并且我们只能够观测到第二轮的抛掷结果，我们用 $y$ 表示第二轮试验的结果，用 $z$ 表示第一轮试验的结果，那么我们称 $y$ 为观测变量，而 $z$ 称为隐变量。最终我们需要根据第二轮的抛掷结果来对参数 $\pi$ ， $p$ 以及 $q$ 进行估计。

对于以上问题EM算法便能够对其进行求解，接下来我们便介绍EM算法的整体思路：

当我们面对一个含有隐变量的概率模型时，我们的目标是极大化观测数据 $Y=(y_1, y_2, ..., y_n)$ 关于参数 $\theta$ 对数似然函数， $Z=(z_1, z_2, ..., z_n)$ 为隐变量，我们需要极大化以下函数：

$L(\theta)=log\;p(Y|\theta)=log\sum_{Z}p(Y, Z|\theta)=log (\sum_{Z}p(Y|Z,\theta)p(Z|\theta) )$

上式运用了全概率公式，而上式就是观测变量的对数似然函数，我们只需对其求极大化便可以求出概率模型中的参数 $\theta$ 。但这一优化问题的求解是困难的。

而求解这一问题，在EM算法中是通过迭代逼近的方式来解决的，假设在第 $i$ 次迭代后 $\theta$ 的估计值为 $\theta^{(i)}$ ，我们所做的事情就是希望在下一次迭代得到的 $\theta$ 能够使得似然函数 $L(\theta)$ 上升，即 $L(\theta)>L(\theta^{(i)})$ ，首先我们计算两者的差值：

$L(\theta)-L(\theta^{(i)})=log (\sum_{Z}p(Y|Z,\theta)p(Z|\theta) )-log\;p(Y|\theta^{(i)})=log(\sum_{Z}p(Z|Y,\theta^{(i)})\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})})-log\;p(Y|\theta^{(i)})$

我们知道对数函数 $log(x)$ 为一个凹函数，则在 $L(\theta)$ 中，我们可以将 $\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})}$ 看做随机变量，而 $p(Z|Y,\theta^{(i)})$ 看做是这个随机变量的概率分布，如下表所示：

$\frac{p(Y\|Z,\theta)p(Z\|\theta)}{p(Z\|Y,\theta^{(i)})}$	$\frac{p(Y\|z_1,\theta)p(z_1\|\theta)}{p(z_1\|Y,\theta^{(i)})}$	$\frac{p(Y\|z_2,\theta)p(z_2\|\theta)}{p(z_2\|Y,\theta^{(i)})}$	...	$\frac{p(Y\|z_n,\theta)p(z_n\|\theta)}{p(z_n\|Y,\theta^{(i)})}$
$p(Z\|Y,\theta^{(i)})$	$p(z_1\|Y,\theta^{(i)})$	$p(z_2\|Y,\theta^{(i)})$	...	$p(z_n\|Y,\theta^{(i)})$

那么

$L(\theta)=log (\sum_{Z}p(Y|Z,\theta)p(Z|\theta) )=log (E(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})}))$

根据Jensen不等式我们知道

$L(\theta)=log (\sum_{Z}p(Y|Z,\theta)p(Z|\theta) )=log (E(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})}))\geq E( log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})}))=\sum_{Z}p(Z|Y, \theta^{(i)})log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})})$

所以：

$L(\theta)-L(\theta^{(i)})\geq \sum_{Z}p(Z|Y, \theta^{(i)})log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})})- log\;p(Y|\theta^{(i)})$

而我们知道：

$\sum_{Z}p(Z|Y, \theta^{(i)})=1$

因此：

$log\;p(Y|\theta^{(i)})=1 * log\;p(Y|\theta^{(i)}) =\sum_{Z}p(Z|Y, \theta^{(i)})log\;p(Y|\theta^{(i)})$

故：

$L(\theta)-L(\theta^{(i)})\geq \sum_{Z}p(Z|Y, \theta^{(i)})log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})})- \sum_{Z}p(Z|Y,\theta^{(i)})log\;p(Y|\theta^{(i)})=\sum_{Z}p(Z|Y, \theta^{(i)})\{log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})})-log\;p(Y|\theta^{(i)})\}=\sum_{Z}p(Z|Y, \theta^{(i)})\{log (\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})p(Y|\theta^{(i)})})\}$

将 $L(\theta^{(i)})$ 移至等号右边得到：

$L(\theta)\geq\sum_{Z}p(Z|Y, \theta^{(i)})\{log(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})p(Y|\theta^{(i)})})\}+L(\theta^{(i)})$

在这里，我们令：

$B(\theta, \theta^{(i)})=\sum_{Z}p(Z|Y, \theta^{(i)})\{log(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})p(Y|\theta^{(i)})})\}+L(\theta^{(i)})$

则

$L(\theta)\geq B(\theta, \theta^{(i)})$

所以 $B(\theta, \theta^{(i)})$ 为观测变量似然函数的下界，我们只需要迭代求取 $B(\theta, \theta^{(i)})$ 的极大值点来近似求取 $L(\theta)$ 的极大值点。迭代公式如下所示：

$\theta^{(i+1)}=arg \;\underset{\theta}{max}B(\theta, \theta^{(i)})$

接着，我们舍去与 $\theta$ 无关的常数部分，得到

$\theta^{(i+1)}=arg \;\underset{\theta}{max}\;B(\theta, \theta^{(i)})=arg \;\underset{\theta}{max}\;\left \{ \sum_{Z}p(Z|Y, \theta^{(i)})\{log(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})p(Y|\theta^{(i)})})\}+L(\theta^{(i)}) \right \}=arg \;\underset{\theta}{max}\;\sum_{Z}p(Z|Y, \theta^{(i)})\{log(\frac{p(Y|Z,\theta)p(Z|\theta)}{p(Z|Y,\theta^{(i)})p(Y|\theta^{(i)})})\}=arg \;\underset{\theta}{max}\;\sum_{Z}p(Z|Y, \theta^{(i)})\{logp(Y|Z,\theta)p(Z|\theta)-logp(Z|Y,\theta^{(i)}p(Y|\theta^{(i)}))\}=arg \;\underset{\theta}{max}\;\left \{\sum_{Z}\{p(Z|Y, \theta^{(i)})logp(Y|Z,\theta)p(Z|\theta)\}-\sum_{Z}\{p(Z|Y, \theta^{(i)})logp(Z|Y,\theta^{(i)}p(Y|\theta^{(i)})\}\right \}=arg \;\underset{\theta}{max}\;\sum_{Z}p(Z|Y, \theta^{(i)})logp(Y|Z,\theta)p(Z|\theta)=arg \;\underset{\theta}{max}\;\sum_{Z}p(Z|Y, \theta^{(i)})logp(Y,Z|\theta)$

我们最后记

$Q(\theta, \theta^{(i)})=\sum_{Z}p(Z|Y, \theta^{(i)})logp(Y,Z|\theta)$

则

$\theta^{(i+1)}=arg \;\underset{\theta}{max}\;Q(\theta, \theta^{(i)})$

实际上我们可以看出， $Q(\theta, \theta^{(i)})$ 函数求的是观测变量与隐变量联合概率分布取对数后关于隐变量的期望值。在每一次迭代的过程中，我们只需要求取使得这个期望值最大的 $\theta$ 即可

综上，EM算法总共两步：

E步：求取 $Q(\theta, \theta^{(i)})$ 函数

M步：求取使得 $Q(\theta, \theta^{(i)})$ 极大的 $\theta$

以下是李航统计学习方法中硬币实验的实现，这里用到之前写过的粒子群算法作为包进行导入，这里对粒子群算法做了很小的一点改动，即run方法中设置返回值，将最优解返回，先给出粒子群算法代码如下(仅供参考)：

import numpy as np
import time
from numpy import random as rd
import matplotlib.pyplot as plt
np.set_printoptions(linewidth=1000, suppress=True)


def objective_func(x, y):
    # 这是一个用来测试的目标函数
    return 2 * np.pi * np.exp(-0.2 * ((x ** 2 + y ** 2) / 2) ** 0.5) + \
            np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) - 2 * np.pi


def objective_func2(x):
    # 这也是一个用来测试的目标函数
    return x ** 2


def GetTime(func_name):
    def wrapper(*args, **kwargs):
        start_time = time.time()
        ret = func_name(*args, **kwargs)
        end_time = time.time()
        run_time = end_time - start_time
        print("pso算法运行时间为%f秒" % run_time)
        return ret
    return wrapper


class PSO(object):

    def __init__(self, alpha, times, obj_type, particle_num, w, c1, c2, obj_func, **kwargs):
        """

        :param alpha: 约束因子，用来控制速度的权重
        :param times: 算法的迭代次数
        :param obj_type: 目标函数的类型，传入"Max"表示目标函数为最大化问题，传入"Min"表示目标函数为最小化问题
        :param particle_num: 粒子数量
        :param w: 惯性因子
        :param c1: 加速常数
        :param c2: 加速常数
        :param obj_func: 被优化的目标函数
        :param kwargs: 目标函数决策变量及取值范围，给定方式例如：x1=(x1_min, x1_max), x2=(x2_min, x2_max), ...,
                        注意此处决策变量及其范围给定的顺序需和目标函数obj_func中决策变量参数顺序一致
        """
        self.alpha = alpha
        self.times = times
        self.obj_type = obj_type
        self.particle_num = particle_num
        self.w = w
        self.c1 = c1
        self.c2 = c2
        self.obj_func = obj_func
        self.kwargs = kwargs
        # self.v_max确定决策变量每一维的最大飞翔速度，根据最大飞翔速度为决策变量取值范围的百分之10到百分之20之间来取值
        self.v_max = (np.array(list(kwargs.values()))[:, 1] - np.array(list(kwargs.values()))[:, 0]) * 0.15
        assert np.all(self.v_max > 0), "输入决策变量的最小值大于最大值"
        self.particles = np.empty((self.particle_num, len(self.kwargs)))
        # for循环初始化初代粒子群
        for dimension in range(len(self.kwargs)):
            var_range = list(self.kwargs.values())[dimension]
            self.particles[:, dimension:dimension + 1] = rd.uniform(var_range[0], var_range[1], \
                                                                    (self.particle_num, 1))
        assert self.obj_type.lower() in ["max", "min"], "目标函数类型填写错误, 请填写'Max'或'Min'"
        if self.obj_type.lower() == "max":
            self.individual_best_fit_value = [-float("inf")] * self.particle_num
        else:
            self.individual_best_fit_value = [float("inf")] * self.particle_num
        self.individual_best_fit_value = np.array(self.individual_best_fit_value)
        self.individual_best_position = np.empty((self.particle_num, len(self.kwargs)))
        # 确定初始速度
        self.v = np.array([rd.uniform(-v_max_value, v_max_value, self.particle_num) \
                           for v_max_value in self.v_max]).T
        # self.bests_pop_fit_values用于存放每次迭代群体最优适应度, self.bests_pop_positions用于存放每次迭代最优群体适
        # 应度对应的位置
        self.bests_pop_fit_values = []
        self.bests_pop_positions = []

    def get_fit_value(self):
        # 求取每个粒子的适应度
        self.fit_value = np.empty(self.particle_num)
        for particle_index in range(self.particle_num):
            decision_var = self.particles[particle_index]
            decision_var_max = np.array(list(self.kwargs.values()))[:, 1]
            decision_var_min = np.array(list(self.kwargs.values()))[:, 0]
            # 一下用于限制决策变量在设定的取值范围之内，如果出了取值范围则设置其适应度为负无穷或正无穷
            if np.all((decision_var_max - decision_var) >= 0) and np.all((decision_var - decision_var_min) > 0):
                self.fit_value[particle_index] = self.obj_func(*tuple(self.particles[particle_index]))
            else:
                if self.obj_type.lower() == "max":
                    self.fit_value[particle_index] = -float("inf")
                else:
                    self.fit_value[particle_index] = float("inf")

    def get_best(self):
        # 获取个体最优以及群体最优粒子位置
        # replace_index为个体最优位置需要替换的粒子的索引
        if self.obj_type.lower() == "max":
            replace_index = self.fit_value > self.individual_best_fit_value
        else:
            replace_index = self.fit_value < self.individual_best_fit_value
        self.individual_best_position[replace_index] = self.particles[replace_index]
        self.individual_best_fit_value[replace_index] = self.fit_value[replace_index]
        # 根据个体最优适应度获得群体最优适应度，以及群体最优粒子位置
        if self.obj_type.lower() == "max":
            self.pop_bset_fit_value = np.max(self.individual_best_fit_value)
        else:
            self.pop_bset_fit_value = np.min(self.individual_best_fit_value)
        best_index = list(self.individual_best_fit_value).index(self.pop_bset_fit_value)
        self.pop_best_position = self.individual_best_position[best_index]
        self.bests_pop_fit_values.append(self.pop_bset_fit_value)
        self.bests_pop_positions.append(self.pop_best_position)

    def calc_v_position(self):
        # 得到粒子群移动后的位置
        self.v = self.v + self.c1 * rd.random() * (self.individual_best_position - self.particles) \
        + self.c2 * rd.random() * (self.pop_best_position - self.particles)
        # 以下两步用于限制速度在[-self.v_max, self.v_max]之内
        self.v = np.array([self.v[r, c] if self.v[r, c] < self.v_max[c] else self.v_max[c] \
                 for r in range(self.particle_num) for c in range(len(self.kwargs))]).reshape(\
                 self.particle_num, len(self.kwargs))
        self.v = np.array([self.v[r, c] if self.v[r, c] > -self.v_max[c] else -self.v_max[c] \
                           for r in range(self.particle_num) for c in range(len(self.kwargs))]).reshape \
        (self.particle_num, len(self.kwargs))
        # 根据计算出的速度获取粒子群移动后的位置
        self.particles = self.particles + self.alpha * self.v

    def run(self):
        for i in range(self.times):
            self.get_fit_value()
            self.get_best()
            self.calc_v_position()
#         print("目标函数最优值为：", self.bests_pop_fit_values[-1])
#         print("最优解为：", self.bests_pop_positions[-1])
        # self.draw()
        return self.bests_pop_positions[-1]

    def draw(self):
        fig = plt.figure()
        ax = plt.subplot(1, 1, 1)
        ax.plot(np.arange(1, self.times + 1), self.bests_pop_fit_values, label="fit_value trend", color="r")
        ax.set_xlabel("iterate times")
        ax.set_ylabel("fit_value")
        ax.grid()
        ax.legend()
        ax.set_title("PSO")
        plt.show()


def main():
    pso = PSO(0.8, 400, "Min", 50, 0.2, 2, 2, objective_func2, x1=(-30, 30))
    pso.run()
    pso = PSO(0.8, 400, "Max", 50, 0.2, 2, 2, objective_func, x1=(-30, 30), x2=(-30, 30))
    pso.run()


if __name__ == "__main__":
    main()

抛掷硬币实验代码如下：

import numpy as np
from numpy import random as rd
from PSO_Algorithm import PSO


class EM(object):
    """
        李航统计学习方法中的模拟硬币实验实现，此处pi, p, q的迭代公式并未采用李航
        统计学习方法中最初给出的公式（统计学习方法第156页）
        而是通过计算Q函数，优化Q函数来求取的，
        此处对于Q函数的优化采用的是之前写过的粒子群算法，直接作为包导入使用的
    """
    
    def __init__(self):
        self.pi_curr, self.p_curr, self.q_curr = [0.4, 0.6, 0.7]
        self.y = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
        print("初始[pi, p, q]:", [self.pi_curr, self.p_curr, self.q_curr])
        print("观测变量：", self.y)
    
    def calc_Q(self, pi, p, q):
        cum_prod_1 = 1
        cum_prod_2 = 1
        for i in range(len(self.y)):
            cum_prod_1 *= (1 - pi) * q ** self.y[i] * (1 - q) ** (1 - self.y[i]) \
                          * (1 - self.pi_curr) * self.q_curr ** self.y[i] * (1 - self.q_curr) \
                          ** (1 - self.y[i]) / \
                          (self.q_curr ** self.y[i] * (1 - self.q_curr) ** (1 - self.y[i]) \
                           * (1 - self.pi_curr) + \
                          self.p_curr ** self.y[i] * (1 - self.p_curr) ** (1 - self.y[i]) * self.pi_curr)
            cum_prod_2 *= pi * p ** self.y[i] * (1 - p) ** (1 - self.y[i]) * self.p_curr ** self.y[i] * \
                          (1 - self.p_curr) ** (1 - self.y[i]) * self.pi_curr / (self.q_curr ** self.y[i] * \
                          (1 - self.q_curr) ** \
                          (1 - self.y[i]) * (1 - self.pi_curr) + self.p_curr \
                          ** self.y[i] * (1 - self.p_curr) ** (1 - self.y[i]) * self.pi_curr)
        return np.log(cum_prod_1) + np.log(cum_prod_2)
    
    def optimize(self):
        pso = PSO(0.8, 250, "Max", 50, 0.2, 2, 2, self.calc_Q, pi=(0, 1), p=(0, 1), q=(0, 1))
        self.pi_curr, self.p_curr, self.q_curr = pso.run()
    
    def run(self):
        for i in range(100):
            self.optimize()
            print("EM算法第%d次迭代[pi, p, q]:" % (i + 1, ), [self.pi_curr, self.p_curr, self.q_curr])

        
def main():
    em = EM()
    em.run()


if __name__ == "__main__":
    main()

$\frac{p(Y\|Z,\theta)p(Z\|\theta)}{p(Z\|Y,\theta^{(i)})}$	$\frac{p(Y\|z_1,\theta)p(z_1\|\theta)}{p(z_1\|Y,\theta^{(i)})}$	$\frac{p(Y\|z_2,\theta)p(z_2\|\theta)}{p(z_2\|Y,\theta^{(i)})}$	...	$\frac{p(Y\|z_n,\theta)p(z_n\|\theta)}{p(z_n\|Y,\theta^{(i)})}$
$p(Z\|Y,\theta^{(i)})$	$p(z_1\|Y,\theta^{(i)})$	$p(z_2\|Y,\theta^{(i)})$	...	$p(z_n\|Y,\theta^{(i)})$