python应用EMD方法 python em算法

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{aligned} 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{aligned}$

其中，$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}] \tag1$

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

1.2 EM算法推导

1.2.1 Jensen 不等式说明

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

若有凸函数$f(x)$,且在函数中取自变量点集$\{x_i\}$,且取对应$\{ \lambda_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)) - log p(y|\theta^i)$

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

$\begin{aligned} 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{aligned}$

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

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

$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)}] \}$

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

$\theta^{(i+1)} = \arg \max_\theta \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) \tag2$

则：

$\theta^{(i+1)} = \arg \max_\theta Q(\theta, \theta^i) \tag3$

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$，且$\mu_j \in \{0,1\},j=1,2,...,n$。即参数$\mu$为模型的隐变量。

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

$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)} \}$

相应的对数似然函数为：

$\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)] \}$

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[log p(y,z|\theta,\theta^{(i)})]$

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

$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)] \}$

1.3.2 M步

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

$\theta^{(i+1)} = \arg \max_\theta 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{aligned} &\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{aligned}$

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

$\begin{aligned} &\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{aligned}$

也即后面两项小于等于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]

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