机器学习之EM算法的原理及推导(三硬币模型)及Python实现

EM算法的简介
EM算法由两步组成：E步和M步，是最常用的迭代算法。

本文主要参考了李航博士的《统计学习方法》
在此基础上主要依据EM算法原理补充了三硬币模型的推导。

1.EM算法的原理

1.1从一个例子开始

三硬币模型
假设有3枚硬币，分别记作A，B和C。 这些硬币正面向上的概率分别是 $\pi ,p$ 和 $q$ 。进行如下抛硬币试验：
1、先抛硬币A, 根据其结果选出硬币B或者硬币C，正面选硬币B，反面选硬币C；
2、然后掷选出的硬币，抛硬币的结果，出现正面记作1，出现反面记作0；
3、独立重复$n$次试验(这里，n=10)，观测结果如下：
$$1,1,0,1,0,0,1,0,1,1$$
假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型参数。

对于单次观测结果，三硬币模型可以写作：
$$\begin{array}{rl}p(y_j|\theta )& =\sum _{z}p(y_j,z|\theta )\\ & =\sum _{z}p(z|\theta )p(y_j|z,\theta )\\ & =\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\end{array}$$

其中，$y_i$是第$j$个观测结果1或0；随机变量$z$是隐变量，表示未观测到的掷硬币A的结果；$\theta =(\pi ,p,q)$是模型参数。
方便起见，观测数据可以表示为$y=(y_1,y_2,...,y_n{)}^T$,隐变量数据可以表示为$z=(z_1,z_2,...,z_n{)}^T$。观测数据的似然函数可以表示为：
$$p(y|\theta )=\sum _{z}p(z|\theta )p(y|z,\theta )$$

即:

$$p(y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]\text{\hspace{1em}(1)}$$

这个问题没有办法直接解析，只能用迭代的方法解决。下面我们先看看EM算法的推导，之后重新再对这个问题进行推导。

1.2 EM算法推导

1.2.1 Jensen 不等式说明

首先说明一下什么是凸函数：
粗糙一点理解，如果函数的二阶导数为正数，那么这个函数就是凸函数：比如开口向上的二次函数就是典型的凸函数。

若有凸函数$f(x)$,且在函数中取自变量点集 { x i } \{x_i\} {xi​},且取对应 { λ i } \{ \lambda_i\} {λi​},满足${\lambda }_i>0,\sum {\lambda }_i=1$,
则有：

$$f(\sum _i{\lambda }_i{x}_i)\le \sum _i{\lambda }_{i}f(x_i)$$

当$x>0$时，$-\log (x)$的二阶导数$\frac{1}{x^2}>0$，故可对$-log(x)$运用Jessen不等式。
如果是对于$log(x)$运用Jessen不等式，不等式方向要变号。

1.2.2 EM算法推导

求解型入（1）式的问题，我们取对数似然函数，也就是对对数似然函数求取极大值：
$$L(\theta )=\log p(y|\theta )=\log (\sum _{z}p(y|z,\theta )p(z|\theta ))$$

运用迭代的思想解决这个问题，假设在第$i$次迭代后$\theta$的估计值是${\theta }^i$。因为想要求取最大对数似然，所以我们希望$L(\theta )>L({\theta }^i)$，并逐步达到极大值，也就是它们的差值达到最大值：

$$L(\theta )-L({\theta }^i)=\log (\sum _{z}p(y|z,\theta )p(z|\theta ))-logp(y|{\theta }^i)$$

利用Jensen不等式，得到其下界:

$$\begin{array}{rl}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)\\ & \ge \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)}-\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)p(y|{\theta }^i)}\end{array}$$

令$J(\theta )=L(\theta )-L({\theta }^i)$,则求取的是：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }J(\theta )$$

J ( θ ) = { ∑ z p ( z ∣ y , θ i ) log ⁡ [ p ( y ∣ z , θ ) p ( z ∣ θ ) ] } − { ∑ z p ( z ∣ y , θ i ) log ⁡ [ p ( z ∣ y , θ i ) p ( y ∣ θ i ) ] } J(\theta) = \{ \sum_z p(z|y,\theta^i) \log [{p(y|z,\theta) p(z|\theta)}] \}- \{ \sum_z p(z|y,\theta^i) \log [{p(z|y,\theta^i) p(y|\theta^i)}] \} J(θ)={z∑​p(z∣y,θi)log[p(y∣z,θ)p(z∣θ)]}−{z∑​p(z∣y,θi)log[p(z∣y,θi)p(y∣θi)]}

因为后一项中无$\theta$项，故：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }\sum _{z}p(z|y,{\theta }^i)\log [p(y|z,\theta )p(z|\theta )]$$

因为：

$$p(y|z,\theta )p(z|\theta )=p(y,z|\theta )$$

设：

$$Q(\theta ,{\theta }^i)=\sum _{z}p(z|y,{\theta }^i)\log p(y,z|\theta )\text{\hspace{1em}(2)}$$

则：

$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^i)\text{\hspace{1em}(3)}$$

EM算法的总结：
E步（求隐变量$p(z|y,{\theta }_i)$）：给定观测数据$y$和当前的参数估计${\theta }_i$,求取隐变量$z$的条件概率分布；
M步：将隐变量当做已知量，求$Q(\theta ,{\theta }_i)$的极大化的$\theta$
E步和M步重复执行，直到收敛。

1.3 三硬币模型继续推导

1.3.1 三硬币模型的隐变量以及完全数据的对数似然函数

我们已知：

$$p(y|\theta )=\prod _{j=1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$

设$y_j$来自掷硬币B的概率为${\mu }_j$, 则来自C的概率为$1-{\mu }_j$，且 μ j ∈ { 0 , 1 } , j = 1 , 2 , . . . , n \mu_j \in \{0,1\},j=1,2,...,n μj​∈{0,1},j=1,2,...,n。即参数$\mu$为模型的隐变量。

于是完全数据的似然函数可以表示为：

p ( y , μ ∣ θ ) = ∏ j = 1 n { [ π p y j ( 1 − p ) 1 − y j ] μ + [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] ( 1 − μ ) } p(y,\mu|\theta) = \prod_{j=1}^n \{ [ \pi p^{y_j} (1-p)^{1-y_j}]^\mu + [(1-\pi)q^{y_j} (1-q)^{1-y_j}]^{(1-\mu)} \} p(y,μ∣θ)=j=1∏n​{[πpyj​(1−p)1−yj​]μ+[(1−π)qyj​(1−q)1−yj​](1−μ)}

相应的对数似然函数为：

log ⁡ p ( y , μ ∣ θ ) = ∑ j = 1 n { μ [ log ⁡ π + y j log ⁡ p + ( 1 − y j ) log ⁡ ( 1 − p ) ] + ( 1 − μ ) [ log ⁡ ( 1 − π ) + y j log ⁡ q + ( 1 − y j ) log ⁡ ( 1 − q ) ] } \log p(y,\mu|\theta) = \sum_{j=1}^n \{\mu [\log \pi + y_j \log p + (1-y_j) \log (1-p)] + (1-\mu) [\log(1-\pi) + y_j \log q + (1-y_j)\log(1-q)] \} logp(y,μ∣θ)=j=1∑n​{μ[logπ+yj​logp+(1−yj​)log(1−p)]+(1−μ)[log(1−π)+yj​logq+(1−yj​)log(1−q)]}

1.3.2 E步：确定$Q$函数

因为EM算法是迭代算法，设第$i$次迭代的参数估计值为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$,又因为隐变量$\mu$代表观测数据来自B的概率，所以第$(i+1)$次隐变量：

$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_i}(1-p^{(i)}{)}^{1-y_i}}{{\pi }^{(i)}(p^{(i)}{)}^{y_i}(1-p^{(i)}{)}^{1-y_i}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_i}(1-q^{(i)}{)}^{1-y_i}}$$

求取$Q$:

$$Q(\theta ,{\theta }_i)=\sum _{z}p(z|y,{\theta }_i)\log p(y,z|\theta )=E_z[logp(y,z|\theta ,{\theta }^{(i)})]$$

将${\mu }_j^{(i+1)}$带入则可以得到：

Q ( θ , θ i ) = ∑ j = 1 n { μ j ( i + 1 ) [ log ⁡ π + y j log ⁡ p + ( 1 − y j ) log ⁡ ( 1 − p ) ] + ( 1 − μ j ( i + 1 ) ) [ log ⁡ ( 1 − π ) + y j log ⁡ q + ( 1 − y j ) log ⁡ ( 1 − q ) ] } Q(\theta, \theta_i)=\sum_{j=1}^n \{\mu_{j}^{(i+1)} [\log \pi + y_j \log p + (1-y_j) \log (1-p)] + (1-\mu_{j}^{(i+1)}) [\log(1-\pi) + y_j \log q + (1-y_j)\log(1-q)] \} Q(θ,θi​)=j=1∑n​{μj(i+1)​[logπ+yj​logp+(1−yj​)log(1−p)]+(1−μj(i+1)​)[log(1−π)+yj​logq+(1−yj​)log(1−q)]}

1.3.2 M步

得到了$Q$函数，接下来就是极大化参数：

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

1.求解$\pi$:

$$\frac{\partial Q(\theta ,{\theta }^i)}{\partial \pi }=\sum _{j=1}^n[{\mu }_j^{(i+1)}\frac{1}{\pi }-(1-{\mu }_j^{(i+1)})\frac{1}{1-\pi }]$$

求取极值，令等式右边为0：

$$\sum _{j=1}^n[{\mu }_j^{(i+1)}\frac{1}{\pi }-(1-{\mu }_j^{(i+1)})\frac{1}{1-\pi }]=0$$

左右两边同时乘$\pi (1-\pi )$得到：

$$\sum _{j=1}^n[{\mu }_j^{(i+1)}(1-\pi )-(1-{\mu }_j^{(i+1)})\pi ]=0$$

$$\sum _{j=1}^n({\mu }_j^{(i+1)}-\pi )=0$$

$$\sum _{j=1}^n{\mu }_j^{(i+1)}-n\pi =0$$

则：
$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$

2.接下来求解$p$:

$$\frac{\partial Q(\theta ,{\theta }^i)}{\partial p}=\sum _{j=1}^n{\mu }_j^{(i+1)}[y_j\frac{1}{p}-(1-y_j^{(i+1)})\frac{1}{1-p}]$$

求取极值，令等式右边为0:

$$\sum _{j=1}^n{\mu }_j^{(i+1)}[y_j\frac{1}{p}-(1-y_j^{(i+1)})\frac{1}{1-p}]=0$$

左右两边同时乘$p(1-p)$得到：

$$\sum _{j=1}^n{\mu }_j^{(i+1)}[y_j(1-p)-(1-y_j^{(i+1)})p]=0$$

$$\sum _{j=1}^n[{\mu }_j^{(i+1)}{y}_j-{\mu }_j^{(i+1)}p]=0$$

则：

$$p^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}$$

3.最后用同样的方法得到$q$:

$$q^{(i+1)}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})}$$

1.3.2 参数空间

1.模型参数
$\pi$: 硬币A正面的概率，在此模型中是一个float类型的数值
$p$: 硬币B正面的概率，在此模型中是一个float类型的数值
$q$:硬币C正面的概率，在此模型中是一个float类型的数值
2.隐变量
$\mu$: 最后观测值到底来源于B还是C，是一个一维向量
$\mu =({\mu }_1,{\mu }_2,...,{\mu }_n)$,其中${\mu }_j$代表第$j$次抛硬币B的概率。

1.4 EM算法的收敛性

证明EM算法的收敛，只需要证明$p(y|{\theta }^{(i)})$是单调递增的即可：

$$p(y|{\theta }^{(i+1)})\ge p(y|{\theta }^{(i)})$$

证明：
由于：

$$p(y|\theta )=\frac{p(y,\theta )}{p(\theta )}\frac{p(y,z,\theta )}{p(y,z,\theta )}=\frac{p(y,z|\theta )}{p(z|y,\theta )}$$

取对数化简得：

$$\begin{array}{rl}& \log p(y|{\theta }^{(i+1)})-\log p(y|{\theta }^{(i)})\\ & =[\log p(y,z|{\theta }^{(i+1)})-\log p(z|y,{\theta }^{(i+1)})]-[\log p(y,z|{\theta }^{(i)})-\log p(z|y,{\theta }^{(i)})]\\ & =[\log p(y,z|{\theta }^{(i+1)})-\log p(y,z|{\theta }^{(i)})]-[\log p(z|y,{\theta }^{(i+1)})-\log p(z|y,{\theta }^{(i)})]\\ & =[\sum _{z}p(z|y,{\theta }^{(i+1)})\log p(y,z|{\theta }^{(i)})-\sum _{z}p(z|y,{\theta }^i)\log p(y,z|{\theta }^{(i)})]-\\ & [\sum _{z}p(z|y,{\theta }^{(i+1)})\log p(z|y,{\theta }^{(i)})-\sum _{z}p(z|y,{\theta }^{(i)})\log p(z|y,{\theta }^{(i)})]\end{array}$$

前两项有$Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})\ge 0$，对后两项进行计算：

$$\begin{array}{rl}& \sum _{z}p(z|y,{\theta }^{(i+1)})\log p(z|y,{\theta }^{(i)})-\sum _{z}p(z|y,{\theta }^{(i)})\log p(z|y,{\theta }^{(i)})\\ & =\sum _z\log [\frac{p(z|y,{\theta }^{(i+1)})}{p(z|y,{\theta }^{(i)})}]p(z|y,{\theta }^{(i)})\\ & \le \log \sum _z\frac{p(z|y,{\theta }^{(i+1)})}{p(z|y,{\theta }^{(i)})}p(z|y,{\theta }^{(i)})\\ & =\log [\sum _{z}p(z|y,{\theta }^{(i+1)})]\\ =0& \end{array}$$

也即后面两项小于等于0，所以$\log p(y|{\theta }^{(i+1)})-\log p(y|{\theta }^{(i)})\ge 0$
得证。

2 三银币模型的Python实现

2.1 模型实现

import numpy as np
np.random.seed(0)


class ThreeCoinsMode(object):
    def __init__(self, n_epoch=5):
        """
        运用EM算法求解三银币模型
        :param n_epoch: 迭代次数
        """
        self.n_epoch = n_epoch
        self.params = {'pi': None, 'p': None, 'q': None, 'mu': None}

    def __init_params(self, n):
        """
        对参数初始化操作
        :param n: 观测样本个数
        :return: 
        """
        self.params = {'pi': np.random.rand(1),
                       'p': np.random.rand(1),
                       'q': np.random.rand(1),
                       'mu': np.random.rand(n)}
        # self.params = {'pi': [0.5],
        #                'p': [0.5],
        #                'q': [0.5],
        #                'mu': np.random.rand(n)}

    def E_step(self, y, n):
        """
        E步：跟新隐变量mu
        :param y: 观测样本
        :param n: 观测样本个数
        :return: 
        """
        pi = self.params['pi'][0]
        p = self.params['p'][0]
        q = self.params['q'][0]
        for i in range(n):
            self.params['mu'][i] = (pi * pow(p, y[i]) * pow(1-p, 1-y[i])) / (pi * pow(p, y[i]) * pow(1-p, 1-y[i]) + (1-pi) * pow(q, y[i]) * pow(1-q, 1-y[i]))

    def M_step(self, y, n):
        """
        M步：跟新模型参数
        :param y: 观测样本
        :param n: 观测样本个数
        :return: 
        """
        mu = self.params['mu']
        self.params['pi'][0] = sum(mu) / n
        self.params['p'][0] = sum([mu[i] * y[i] for i in range(n)]) / sum(mu)
        self.params['q'][0] = sum([(1-mu[i]) * y[i] for i in range(n)]) / sum([1-mu_i for mu_i in mu])

    def fit(self, y):
        """
        模型入口
        :param y: 观测样本
        :return: 
        """
        n = len(y)
        self.__init_params(n)
        print(0, self.params['pi'], self.params['p'], self.params['q'])
        for i in range(self.n_epoch):
            self.E_step(y, n)
            self.M_step(y, n)
            print(i+1, self.params['pi'], self.params['p'], self.params['q'])

2.2 模型测试结果

def run_three_coins_model():
    y = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
    tcm = ThreeCoinsMode()
    tcm.fit(y)

运行结果：

0 [0.5488135] [0.71518937] [0.60276338]
1 [0.54076424] [0.65541668] [0.53474516]
2 [0.54076424] [0.65541668] [0.53474516]
3 [0.54076424] [0.65541668] [0.53474516]
4 [0.54076424] [0.65541668] [0.53474516]
5 [0.54076424] [0.65541668] [0.53474516]

参考文献：
《统计学习方法》 李航著