# EM算法与GMM模型 -- 潘登同学的Machine Learning笔记

GMM模型
- 单高斯模型 GM的参数估计（本质是最大似然估计）
- 混合高斯分布 GMM 的参数估计
EM算法
GMM算法应用在图像识别中

GMM模型

高斯密度函数估计是一种参数化模型。高斯混合模型（Gaussian Mixture Model, GMM）是单一高斯概率密度函数的延伸，GMM能够平滑地近似任意形状的密度分布。高斯混合模型种类有单高斯模型（Single Gaussian Model, SGM）和高斯混合模型（Gaussian Mixture Model, GMM）两类。类似于聚类，根据高斯概率密度函数（Probability Density Function, PDF）参数不同，每一个高斯模型可以看作一种类别，输入一个样本x，即可通过PDF计算其值，然后通过一个阈值来判断该样本是否属于高斯模型。很明显，SGM适合于仅有两类别问题的划分，而GMM由于具有多个模型，划分更为精细，适用于多类别的划分，可以应用于复杂对象建模。

单高斯模型 GM的参数估计（本质是最大似然估计）

具体用法是这样，先用训练样本估计出单高斯模型的参数后，再用这个高斯模型去判断一个样本是否属于该类别。

输入: 给定一组样本 $X_1,X_2,\ldots,X_n$ ，已知它们来自于高斯分布 $N(\mu,\sigma)$ 。
输出： $\mu,\sigma$

单个样本 $X_k$ 的似然函数（表示 $X_k$ 出现的概率）
$f(X_k) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(X_k-\mu)^2}{2\sigma^2}}$

$X_1,X_2,\ldots,X_n$ 全部发生的概率
$L(x) = \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(X_i-\mu)^2}{2\sigma^2}}$

取对数，变连乘为连加
$\begin{aligned} l(X) &= \sum_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(X_i-\mu)^2}{2\sigma^2}}\\ &= -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n(X_i-\mu)^2\\ \end{aligned}$

分别对 $\mu 和 \sigma$ 求偏导，令偏导数=0
$\mu = \frac{1}{n}\sum_{i=1}^n X_i\\ \sigma^2 = \frac{1}{n}\sum_{i=1}^n(X_i-\mu)^2\\$

最后的结果，其实就是一个样本均值和样本方差；(从某种程度上也说明最大似然估计出来的结果是正确的)

混合高斯分布 GMM 的参数估计

举个栗子

考虑这样一个问题, xx市的教育局对两所学校（A学校 & B学校）抽查期末概率论的成绩，但是粗心的数据分析师把潘登同学的学校那一栏删掉了（为什么受伤的是我）。现在根据已有的数据能不能找出潘登同学是哪个学校的呢？

回答：很简单，对已有的数据（A学校、B学校）分别进行单高斯模型的参数估计，然后再把潘登同学的分数带进模型中，求似然，取大的那个就是正确答案了

进一步地，如果这个粗心的数据分析把所有同学的学校那一栏都删掉了，但是他还想知道潘登同学属于哪个学校的？（已知A学校的均分高于B学校）

这个答案是肯定的，但是在解决之前要对数据做一项假设：

这些数据来源于两个不同的正态总体

假如所有学生的成绩服从如下分布:

混合高斯模型

那么我们只要能把两个正态总体区分开来，估计出不同正态分布的参数，再求似然不就能得出结果了嘛...

GMM 混合高斯分布

假设随机变量 X 是由 K 个高斯分布混合而来，取到各个高斯分布的概率为 $\Pi_1, \Pi_2 ,\ldots, \Pi_K$ ，第 i 个高斯分布的均值为 $μ_i$ ，标准差为 $σ_i$ ，根据只观测到一系列样本 $X_1,X_2,\ldots,X_n$ ，估计 $\Pi,\mu,\sigma$

注意 这里的 $\Pi$ 应该理解为先验概率，根据前面的例子就是，每个学校分别都有多少人，假如说A学校人数：B学校人数 = 2 ：8，那么 $\Pi_A = 0.2, \Pi_B = 0.8$

GMM 模型的似然函数

根据前面同样的操作，得出 GMM 模型的似然函数
$l_{\pi,\mu,\sigma} = \sum_{i=1}^n \log(\sum_{k=1}^K\pi_kN(x_i|\mu_k,\sigma_k))$

分两步求解 GMM

第0步：随机初始化 $\Pi,\mu,\sigma$
第一步：估计数据来源于哪个分布
$\gamma(i,k) = \frac{\pi_k N(x_i|\mu_k,\sigma_k)}{\sum_{j=1}^K\pi_jN(x_i|\mu_j,\sigma_j)}$

解释一下：这个 $\gamma(i,k)$ 表示的是 $P(X\in Z_k|X=X_i)$ 是一种条件概率，其实他是由贝叶斯公式来的；

我们知道贝叶斯公式沟通了先验概率与后验概率:
$P(A|B) = \frac{P(B|A)P(A)}{P(B)}$

而 $P(A)$ 可以理解为 $\pi_k$ 表示这一类别出现的概率, 而 $P(B|A)$ 可以理解为 $N(x_i|\mu_k,\sigma_k)$ 表示在这一类别出现的条件下 $x_i$ 出现的概率；而 $P(B)$ 则表示 $x_i$ 出现的全概率，即 $\sum_{j=1}^K\pi_jN(x_i|\mu_j,\sigma_j)$ 可以理解为 $x_i$ 的全概率公式；
而 $P(A|B)$ 就表示为， $x_i$ 属于这一类别的概率，显然 $\sum_{i=1}^nP(X\in Z_k|X=X_i) = 1$

这样的跟之前讲的softmax多分类很像，那我们只要固定住 $i$ ,求出最大 $\gamma(i,k)$ ，那么 $i$ 就属于第 $k$ 类。

第二步：更新 $\Pi,\mu,\sigma$

根据第一步求解出的 $\gamma(i,k)$ ,来将样本 $x_i$ 划分到k个高斯分布中。更新 $\Pi,\mu,\sigma$
$\Pi_k = \sum_{i=1}^n\gamma(i,k)\\ \mu_k = \frac{1}{\Pi_k}\sum_{i=1}^n\gamma(i,k)x_i\\ \sigma_k = \frac{1}{\Pi_k}\sum_{i=1}^n\gamma(i,k)(x_i-\mu_k)^2$

算法总结

0.随机初始化 $\Pi,\mu,\sigma$
1.计算 $\gamma(i,k)$
2.根据 $\gamma(i,k)$ 更新 $\Pi,\mu,\sigma$
3.重复1、2，直到 $\Pi,\mu,\sigma$ 不再变化

回到前面的栗子

我们假设A学校的均分是高于B学校的，因为总要有一个均分高，不妨设为A，根据成绩分布的概率密度函数，模拟出所有学生的成绩数据，再根据GMM模型，算出潘登同学应该属于哪个学校。

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
from sklearn.mixture import GaussianMixture
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def pdf_normal(mu_1, sigma_1, mu_2, sigma_2, plot=True):
    '''
    mu_1: mu_1小于mu_2.
    生成混合高斯分布的概率密度曲线.
    '''
    x = np.linspace(mu_1-3*sigma_1, mu_2+3*sigma_2, 100)
    pro_1 = st.norm.pdf(x, mu_1, sigma_1)
    pro_2 = st.norm.pdf(x, mu_2, sigma_2)
    pro = (pro_1 + pro_2)/(pro_1 + pro_2).sum()
    if plot:
        plt.plot(x, pro, 'r-')
        plt.ylabel('概率密度')
        plt.title(f'均值为{mu_1}标准差为{sigma_1} 与 均值为{mu_2}标准差为{sigma_2}的和的概率密度曲线')
        plt.show()
    return x,pro
    
# def GMM(x,)
    
if __name__=='__main__':
    np.random.seed(42)
    x,pro = pdf_normal(90, 1.4, 95, 1, False)
    grade = []
    for i,j in zip(x,pro):
        for k in range(round(j*1e5)):
            grade.append(i)
    gmm = GaussianMixture(n_components=2, covariance_type="tied")
    grade = np.array(grade).reshape(-1,1)
    gmm = gmm.fit(grade)
    pd_grad = np.array([[96]])
    cluster = gmm.predict(pd_grad)
    print('两个学校的均分分别是：\n',gmm.means_)
    school = {0:'A',1:'B'}
    print('潘登同学属于学校：', school[int(cluster)])
    pred = gmm.predict(grade)
    plt.hist(grade[pred==0], color='red', density=True,bins=100,label='B学校的成绩')
    plt.hist(grade[pred==1], color='blue', density=True,bins=100,label='A学校的成绩')
    plt.legend()
    plt.show()

结果如下：

潘登同学属于哪个学校

成绩归属

我们可以看到，最终GMM模型得出的均值，与我们模拟数据的两个正态总体的均值90与95非常接近，然后潘登同学96分，来自A学校的概率更大些，所以最终结果也是A学校。

EM算法

EM算法是一种迭代优化策略，由于它的计算方法中每一次迭代都分两步，其中一个为期望步（E步），另一个为极大步（M步），所以算法被称为EM算法（Expectation-Maximization Algorithm）。EM算法受到缺失思想影响，最初是为了解决数据缺失情况下的参数估计问题，

其基本思想是：首先根据己经给出的观测数据，估计出模型参数的值；然后再依据上一步估计出的参数值估计缺失数据的值，再根据估计出的缺失数据加上之前己经观测到的数据重新再对参数值进行估计，然后反复迭代，直至最后收敛，迭代结束。

EM算法的步骤

由单一分布的似然估计推出多分布的似然估计

对于n个样本数据 $x = (x_1,x_2,\ldots,x_n)$ , 模型参数为 $\theta$ , 最大化模型的对数似然
$\hat{\theta} = \argmax\sum_{i=1}^n\log p(x_i;\theta)$

如果我们观测到数据有未观察到的隐含数据 $Z=(Z_1,Z_2,\ldots,Z_k)$ ,即上面的n个样本归属于哪一类是未知的，则改写为：
$\begin{aligned} \hat{\theta} &= \argmax\sum_{i=1}^n\log p(x_i;\theta)\\ &= \argmax\sum_{i=1}^n\log\sum_{j=1}^kp(x_i,z_j;\theta)\\ &= \argmax\sum_{i=1}^n\log\sum_{j=1}^k Q_i(z_j)\frac{p(x_i,z_j;\theta)}{Q_i(z_j)}\\ &\geq \argmax\sum_{i=1}^n\sum_{j=1}^kQ_i(z_j)\log\frac{p(x_i,z_j;\theta)}{Q_i(z_j)} (*)\\ \end{aligned}$

注意:

这里引入了一个新的分布 $Q_i(z_j)$
$(*)$ 式是由Jensen不等式得到的：

如果f是凸函数，X是随机变量，那么: $E[f(X)]\geq f(E(X))$ 。当且仅当X是常量时，该式取等号。其中，E(X)表示X的数学期望。

例如： $f(\theta x + (1-\theta)y)\leq\theta f(x) + (1-\theta)f(y)$

若 $\theta_1,\theta_2,\ldots,\theta_n \geq 0,\theta_1 + \cdots + \theta_k = 1$ , 则 $f(\theta_1 x_1 + \cdots + \theta_n x_n)\leq\theta_1 f(x_1) + \cdots + \theta_n f(x_n)$
由多分布的似然估计求后验概率

而 $(*)$ 式中log是凸函数，把 $\frac{p(x_i,z_j;\theta)}{Q_i(z_j)}$ 理解为随机变量，当Jensen不等式取等号时，有
$\frac{p(x_i,z_j;\theta)}{Q_i(z_j)} = c,c为常数$

由于 $Q_i(z_j)$ 是一个分布，所以满足 $\sum_{j=1}^k Q_i(z_j) = 1$ ，则 $\sum_{j=1}^k p(x_i,z_j;\theta) = c$ ,则有
$Q_i(z_j) = \frac{p(x_i,z_j;\theta)}{\sum_{j=1}^k p(x_i,z_j;\theta)} = \frac{p(x_i,z_j;\theta)}{p(x_i;\theta)} = p(z_j|x_i;\theta)$

这里也是用到了全概率公式、贝叶斯公式，来计算后验概率。

(算法)EM算法

0.随机初始化模型参数 $\theta$
1.E步：计算联合分布的条件概率期望
$Q_i(z_j) = p(z_j|x_i;\theta_t)\\ l(\theta_t) = \sum_{i=1}^n\sum_{j=1}^kQ_i(z_j)\log\frac{p(x_i,z_j;\theta_t)}{Q_i(z_j)}$
2.M步：极大化 $l(\theta_t)$ ,得到 $\theta_{t+1}$
$\theta_{t+1} = \argmax l(\theta_t)$
3.重复23步直到 $\theta_t$ 已经收敛，算法结束

应用EM框架

把EM框架用在GMM上

现在对比EM算法于GMM模型，EM算法中的 $Q_i(z_j)$ 其实就是GMM的 $\gamma(i,k)$ ；

GMM的三步也是先随机，再计算条件概率，再更新参数；所以说GMM的本质是EM用于估计混合高斯分布的情况；

把EM框架用在K-means上

回想K-means，也是三步，先随机中心点，再计算其他点到中心点距离，归堆更新中心点；其实也可以看作使用了EM算法的框架；

GMM算法应用在图像识别中

GMM基于的假设是不同类别的东西应该服从不同的正态分布，也就是说一张图片的各个部位应该应该服从不同的正态分布，依据不同的正态分布将图片的各个部分分开；

观测如下一朵花,我们想把花瓣与背景进行分离。

Flower

应用GMM模型：

from sklearn.mixture import GaussianMixture
from sklearn.cluster import AgglomerativeClustering
import matplotlib.pyplot as plt
from PIL import Image
import numpy as np

np.random.seed(423)
im = Image.open('flower1.jpg')
plt.imshow(im)
newdata = np.array(im).reshape((-1, 3))
shape = newdata.shape
print(shape)
gmm = GaussianMixture(n_components=3, covariance_type="tied")
#高斯分布主要是依据颜色的分布进行划分，n_components越小，划分的越粗糙
gmm = gmm.fit(newdata)

cluster = gmm.predict(newdata)
# cluster = gmm.predict_proba(newdata)
# cluster = (cluster > 0.98).argmax(axis=1)
cluster = cluster.reshape(842, 474)
print(cluster.shape)
plt.imshow(cluster)
plt.show()

结果如下:

Flower_afterGMM

EM算法与GMM模型就是这样了, 继续下一章吧！pd的Machine Learning