EM算法和Python代码实现

reept

已于 2024-04-10 17:45:39 修改

阅读量1k

点赞数 16

文章标签：概率论机器学习人工智能 python

于 2024-02-01 17:10:07 首次发布

本文链接：https://blog.csdn.net/reept/article/details/135819177

版权

1、简介

EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。

一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z都具备，则称为完全数据，观测数据Y又称为不完全数据，假设给定观测数据Y，其概率分布是P（Y|θ），其中θ是要估计的模型参数，那么不完全数据Y的似然函数是logP(Y|θ)，假设Y和Z的联合概率分布是P（Y，Z|θ）,则完全数据的对数似然函数是logP(Y,Z|θ）。

EM算法由极大似然估计推导得出，一般过程就是从似然函数L(θ)出发进行如下推导：

$L(\theta ) = logP(Y\mid \theta ) = log \sum_{Z}^{} P(Y,Z \mid \theta ) = Log(\sum_{Z}^{}P(Y\mid Z,\theta )P(Z\mid \theta ))$

然后通过Jesen不等式得到其下界 $B(\theta ,\theta ^{(i)})$ ：

通过求解使得 $B(\theta ,\theta ^{(i)})$ 最大的 $\theta ^{i+1}$ 进行迭代，求解 $argmax_{\theta }B(\theta ,\theta ^{i})$ = = 因此，求 $B(\theta ,\theta ^{(i)})$ 的最大化参数θ换成了求Q函数的最大化参数。

Q函数是完全数据的对数似然函数 $logP(Y,Z\mid \theta )$ 关于在给定观测数据Y和当前参数 $\theta ^{i}$ 下对未观测数据Z的条件概率分布 $P(Z\mid Y,\theta ^{i})$ 的期望,即：

EM算法通过初设参数 $\theta ^{i}$ ，求出Q函数，作为E步，然后通过求解在E步Q函数基础上的最大化参数 $\theta ^{i+1}$ ,作为M步，然后重复E、M步直到收敛为止。

EM算法：

输入：观测变量数据Y，隐变量数据Z（这里数据并不是Z的观测数据，而是Z的可能取值），联合分布P(Y,Z|θ），条件分布P（Z|Y,θ）；

输出：模型参数θ；

(1）选择参数的初值 $\theta ^{(0)}$ ,开始迭代；

（2）E步：记 $\theta ^{(i)}$ 为第i次迭代参数θ的估计值，在第i+1次迭代的E步，计算：

这里 $P(Z|Y,\theta ^{(i)})$ 是在给定观测数据Y和当前的参数估计参数估计 $\theta ^{(i)}$ 下隐变量数据Z的条件概率分布；

（3）M步求使得 $Q(\theta ,\theta ^{(i)})$ 极大化的θ，确定第i+1 次迭代的参数的估计值 $\theta ^{(i+1)}$

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

(4) 重复第（2）步和第（3）步，直到收敛。

关于EM算法的几点说明：

1、参数的处置可以任意选择，但EM算法对初值是敏感的，因此选择初值十分重要，常用的办法是选择几个不同的初值进行迭代，然后对得到的各个参数估计值进行评估，从中择优；

2、E步求 $Q(\theta ,\theta ^{(i)})$ ，始终Z是未观测数据，Y是观测数据，第一个参数θ表示本轮需要极大化的参数，而第二个参数表示上一轮已经求得的参数或者 $\theta ^{(0)}$ ，每次迭代实际在求Q函数及其极大。

3、M步求 $Q(\theta ,\theta ^{(i)})$ 极大化，得到 $\theta ^{(i+1)}$ ,完成 $\theta ^{(i)}$ -> $\theta ^{(i+1)}$ 。每次求 $Q(\theta ,\theta ^{(i)})$ ，实际就是使似然函数增大或达到局部极值的过程。

4、给出停止迭代的条件，一般是对较小 $\varepsilon _{1}$ 、 $\varepsilon _{2}$ ，满足待：

$\left \| \theta ^{i+1} - \theta ^{i}\right \| < \varepsilon _{1}$ 或者 $\left \| Q(\theta ^{i+1 },\theta ^{i}) - Q(\theta ^{i },\theta ^{i}) \right \| < \varepsilon _{2}$ 则停止迭代。

EM算法的应用范围很广，一个重要的应用是高斯混合模型的参数估计，是高斯混合模型参数估计的有效方法。

2、高斯混合模型及其参数估计的EM算法

高斯混合模型是指有如下形式概率密度函数的模型：

模型生成数据的过程：首先依概率 $\alpha _{k}$ 选择第k个模型 $\Phi(y\mid \theta _{k} )$ ,然后依据 $\Phi(y\mid \theta _{k} )$ 生成观测数据 $y_{j}$ j= 1,2, 3...N。反映观测数据 $y_{j}$ 来自第k个模型的数据 $\gamma _{jk}$ 是未知的。 $\gamma _{jk}$ 的定义如下：

$\gamma _{jk}$ 是0-1随机变量。

有了以上基本定义，可以看出模型的参数θ = { $\alpha$ 、 $\mu$ 、 $\sigma ^{2}$ } ，（ $\mu$ 表示分模型的均值， $\sigma ^{2}$ 表示分模型的方差， $\alpha$ 表示分模型的命中概率，都是长度为k的向量）隐变量： $\gamma _{jk}$ 。就可以推导出高斯混合模型的EM算法。

1）、首先求出完全数据的似然函数 $P(Y,\gamma \mid \theta )$ ,最后求出完全数据的对数似然函数：

2）、然后进行EM算法的E步，确定Q函数：

$Q(\theta ,\theta ^{(i)})$ = $E_{\gamma }\left [ logP(y,\gamma \mid \theta )\mid y,\theta ^{i+1} \right ]$ = $E_{\gamma }\left [P(\gamma \mid Y,\theta ^{i}) logP(Y,\gamma \mid \theta ) \right ]$

由于 $P(\gamma \mid Y,\theta ^{i} )$ = $\frac{P(\gamma ,Y \mid\theta ^{i})}{P(Y|\theta^{i} )}$ 代入上式可得 $Q(\theta ,\theta ^{(i)})$ = $E_{\gamma }\left [\frac{P(\gamma ,Y \mid\theta ^{i})}{P(Y|\theta^{i} )} logP(Y,\gamma \mid \theta ) \right ]$

由于 $\theta ^{i}$ 已知，所以上式左边分子式为已知，从而：

$E_{\gamma }\left [\frac{P(\gamma ,Y \mid\theta ^{i})}{P(Y|\theta^{i} )} logP(Y,\gamma \mid \theta ) \right ] \Leftrightarrow E_{\gamma }\left [logP(Y,\gamma \mid \theta ) \right ]$

因而 $Q(\theta ,\theta ^{(i)})$ 可以去上式右边的式子，即

$Q(\theta ,\theta ^{(i)})$ = $E_{\gamma }\left [logP(Y,\gamma \mid \theta ) \right ]$

上式就是对log $P(Y,\gamma \mid \theta )$ ）求已知γ条件概率分布的情况下的期望（将θ参数视为常量）。因此得到：

$Q(\theta ,\theta ^{(i)})$

其中 $E(\gamma _{jk})$ 记为 $\gamma \hat{}_{jk}$

表示当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据 $y_{j}$ 的响应度。

3）、确定EM算法的M步

迭代的M步是求函数 $Q(\theta ,\theta ^{(i)})$ 对θ的极大值，

重复以上计算，直到对数似然函数值不再有明显的变化为止。

3、Python实现算法

实现算法可参考以下链接：

1、numpy实现周志华机器学习 9.4.3 高斯混合聚类（GMM算法）_多维正态分布-CSDN博客

2、机器学习之高斯混合模型(GMM)及python实现_gmm python-CSDN博客

（以上实现的是多维高斯分布混合模型）

对于第二链接，模型fit时候，只有迭代次数，没有给出收敛条件（博客留出的问题），fit()函数内收敛条件代码：

 def fit(self, y, max_iter=100):
        y = np.array(y)
        self.N, self.D = y.shape
        self.__init_params()
        Epsilon = 0.001
        count = 0
        for _ in range(max_iter):
            Old_Mu,Old_Sigma,Old_alpha = copy.deepcopy(self.params['mu']),copy.deepcopy(self.params['Sigma']),copy.deepcopy(self.params['alpha'])
            
            Old_Mu,Old_Sigma,Old_alpha = np.sum(self.params['mu'],axis = 0),np.sum(self.params['Sigma'],axis = 0),Old_alpha
            self._E_step(y)
            self._M_step(y)
            print(self.params['mu'].shape)
            count  +=1
            
            if abs(np.linalg.norm(x=np.sum(self.params['mu'],axis = 0), ord=2) - np.linalg.norm(x=Old_Mu, ord=2)) < Epsilon and abs(np.linalg.norm(x=np.sum(self.params['Sigma'],axis = 0), ord=2) - np.linalg.norm(x=Old_Sigma, ord=2))  < Epsilon and abs(np.linalg.norm(x=self.params['alpha'], ord=2) - np.linalg.norm(x=Old_alpha, ord=2)) < Epsilon:  
                print(count)
                break

代码中用到的条件是参数θ收敛条件，即 $\left \| \theta ^{i+1} - \theta ^{i}\right \| < \varepsilon _{1}$ ，这里取得是前后两次参数（代码中用）按特征维度求和，得到一个和向量，再求这个和向量2范数的差值绝对值。

代码运行效果：

图1 代码运行效果

图2 代码运行效果

其中55 是迭代停止的次序数，原代码中max_iter =100, 矩阵表示 $\gamma _{jk}$ =1 时，数据所属的高斯分布。聚类图表示模型适配完成时或参数更新完成时，二维数据所属高斯分布（用不同颜色表示不同高斯分布，或者类别，分成三类）。

4、参考资料

（1）《机器学习》周志华

（2）《统计学习》李航

（3）numpy实现周志华机器学习 9.4.3 高斯混合聚类（GMM算法）_多维正态分布-CSDN博客

（4）机器学习之高斯混合模型(GMM)及python实现_gmm python-CSDN博客

reept

关注

16
点赞
踩
20

收藏

觉得还不错? 一键收藏
1
评论
EM算法和Python代码实现

EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z都具备，则称为完全数据，观测数据Y又称为不完全数据，假设给定观测数据Y，其概率分布是P（Y|θ），其中θ是要估计的模型参数，那么不完全数据Y的似然函数是logP(Y|θ)，假设Y和Z的联合概率分布是P（Y，Z|θ）,则完全数据的对数似然函数是logP(Y,Z|θ）。
复制链接

扫一扫