经典机器学习模型(九)EM算法在高斯混合模型中的应用

EM算法的推导可以参考：

经典机器学习模型(九)EM算法的推导

若随机变量X服从一个数学期望为$\mu$、方差为${\sigma }^2$的正态分布，可以记为$N(\mu ，\sigma 2)$。

如果我们知道一组样本服从高斯分布，那么根据极大似然估计，就很容得出参数$u,{\sigma }^2$的值，求解过程如下：

注：

上图用均值而不是平均值，这两者的差别在于：期望还可以称为均值，是对应于总体分布的，但是平均值一般指用样本计算的算术平均数。

高斯混合模型，既然是高斯混合，那么一定是要有高斯分布的，接着还有另一个分布和它混合在一起。

假如我们有一个标准正态分布，即均值为0，方差1，另一个高斯分布的均值为5，方差为1，我们按照1:1叠加合并，也就是分别对应的占比是0.5，就得到了下图所示的混合高斯分布。

假如，我们现在有K个高斯分布：
$$\begin{gathered} N(u_1,{\sigma }_1^2),N(u_2,{\sigma }_2^2),...,N(u_K,{\sigma }_K^2) \\ 对应的比例为：a_1,a_2,...,a_K,意味着我们给每个高斯分布设一个权重 \\ 比例之和为1，即a_1+a_2+...+a_K=1 \end{gathered}$$
如果我们想要估计出混合高斯分布里面每一个分布的参数$u$和${\sigma }^2$，此时就麻烦了。因为我们只能观测到样本，并不知道哪一个样本属于哪一个高斯分布。

如果，我们知道每一个样本属于哪一个高斯分布，那么就很简单了，利用极大似然估计，就能直接求出每一个高斯分布所拥有的样本的均值和方差，因此关键点就在于找出每一个样本属于哪一个高斯分布。

1 高斯混合模型公式推导

1.1 高斯混合模型E步推导

根据上文，我们现在并不知道对应的样本到底是哪个高斯分布，这就是对应的隐变量。

这恰好就是EM算法能够实现的情况。

所以，我们就先准备高斯混合模型的E步，这一步要完成的是两个目的：

• 一是补全信息，也就是下文中即将要出场的隐变量${\gamma }_{jk}$；
• 二是确定 M 步的目标函数。

1.1.1 隐变量的表示

$$\begin{gathered} 我们设{\gamma }_{jk}=\{\begin{array}{ll}1,& y_j来自第k个高斯分布N(u_k,{\sigma }_k^2)\\ 0,& \text{其他}\end{array} \\ 例如观测的样本来自于第1个高斯分布，可以表示为： \\ {\gamma }_{j1}=\{\begin{array}{ll}1,& 观测样本y_j来自第1个高斯分布N(u_1,{\sigma }_1^2)\\ 0,& \text{其他}\end{array} \\ 那么，此时我们有y_1,y_2,...,y_N观测样本，那么隐变量就可以表示为： \\ {y}_1,{\gamma }_{11},{\gamma }_{12},{\gamma }_{13},...,{\gamma }_{1K} \\ {y}_2,{\gamma }_{21},{\gamma }_{22},{\gamma }_{23},...,{\gamma }_{2K} \\ ...... \\ {y}_N,{\gamma }_{N1},{\gamma }_{N2},{\gamma }_{N3},...,{\gamma }_{NK} \\ 对应的y_1,y_2,...,y_N就是「观测量」，是已知的。而后面的「隐变量」{\gamma }_{ji}自然是未知的。 \\ 此时对应的参数\theta = \\ [a_1,u_1,{\sigma }_1^2 \\ {a}_2,u_2,{\sigma }_2^2 \\ ...... \\ {a}_K,u_K,{\sigma }_K^2] \end{gathered}$$

1.1.2 隐变量的求解

${\gamma }_{jk}$是在当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，那么如何求出${\gamma }_{jk}$呢？

我们可以写成期望的形式：
$$\begin{gathered} {\gamma }_{jk}=E({\gamma }_{jk}|y_j,{\theta }^i) \\ =1\ast P({\gamma }_{jk}=1|y_j,{\theta }^i)+0\ast P({\gamma }_{jk}=0|y_j,{\theta }^i) \\ =P({\gamma }_{jk}=1|y_j,{\theta }^i) \end{gathered}$$
然后，我们借助贝叶斯公式：
$$\begin{gathered} {\gamma }_{jk}=P({\gamma }_{jk}=1|y_j,{\theta }^i) \\ =\frac{P(y_j|{\gamma }_{jk}=1,{\theta }^i)P({\gamma }_{jk}=1|{\theta }^i)}{P(y_j,{\theta }^i)} \\ 再利用全概率公式p(A)=\sum p(A|B_i)p(B_i)，将分母展开 \\ =\frac{P(y_j|{\gamma }_{jk}=1,{\theta }^i)P({\gamma }_{jk}=1|{\theta }^i)}{{\sum }_{k=1}^{K}P(y_j|{\gamma }_{jk}=1,{\theta }^i)P({\gamma }_{jk}=1|{\theta }^i)} \\ 上式中{\theta }^i是上一轮的参数，是已知的。 \\ 观察上式，P(y_j|{\gamma }_{jk}=1,{\theta }^i)是指高斯分布参数{\theta }^i已知， \\ 样本y_j属于第k个分布，那么可以用高斯密度函数进行表示\varnothing (y_j|u_k^i,{\sigma }_k^{2(i)}) \\ P({\gamma }_{jk}=1|{\theta }^i)表示已知上一轮的参数{\theta }^i，那么属于第k个分布的概率很容易得出为a_k^{(i)} \\ 因此，我们可以求得隐变量的值，如下式： \\ {\gamma }_{jk}=\frac{\varnothing (y_j|u_k^i,{\sigma }_k^{2(i)})a_k^{(i)}}{{\sum }_{k=1}^K\varnothing (y_j|u_k^i,{\sigma }_k^{2(i)})a_k^{(i)}} \\ 我们令n_k=\sum _{j=1}^N{\gamma }_{jk},表示对所有的样本进行统计，如果属于第k个分布，那么{\gamma }_{jk}=1 \\ {n}_k的意思就是属于第k个分布的样本数 \end{gathered}$$
根据以上推导的公式，我们就能在已知上一轮参数${\theta }^{(i)}$的情况下，求得隐变量${\gamma }_{jk}$的值。

注：

我们现在在已知上一轮参数${\theta }^{(i)}$的情况下，求出了隐变量${\gamma }_{jk}$的值，其实我们现在已经知道哪一个样本属于哪一个高斯分布了，因此只要求得每一个高斯分布所有样本的均值和方差，就求出了下一轮的的参数${\theta }^{(i+1)}$。
$$\begin{gathered} 我们现在求出了隐变量{\gamma }_{jk}的值,属于第k个高斯分布的样本数我们表示为n_k=\sum _{j=1}^N{\gamma }_{jk} \\ 下一轮的a_k是第k个高斯分布样本所占的比例，很容易得出：a_k=\frac{n_k}{N} \\ 下一轮的u_k是第k个高斯分布的均值,很容易得出：u_k=\frac{{\sum }_{j=1}^N{r}_{jk}{y}_j}{n_k} \\ 下一轮的{\sigma }^2是第k个高斯分布的方差,很容易得出：{\sigma }_k^2=\frac{{\sum }_{j=1}^N{r}_{jk}(y_j-u_k{)}^2}{n_k} \end{gathered}$$
我们下面通过另一个角度，即通过联合概率分布构造目标函数$Q(\theta ,{\theta }^{(i)})$进行求解。

1.1.3 联合概率分布

此时的联合概率分布，因为是连续数据，所以我们可以用联合概率密度来表示，对应的值应该是连乘的形式：

然接着我们生成完全数据的对数似然函数：
$$\begin{gathered} lnL(\theta )=\sum _{k=1}^K(n_{k}ln{a}_k+\sum _{j=1}^{N}ln[\varnothing (y_j|{\theta }_k){]}^{r_{jk}}) \\ =\sum _{k=1}^K(n_{k}ln{a}_k+\sum _{j=1}^N{r}_{jk}ln[\varnothing (y_j|{\theta }_k)]) \\ 其中n_k=\sum _{j=1}^N{\gamma }_{jk} \end{gathered}$$
根据以上推导，我们从E步就得到了目标函数为：
$$\begin{gathered} Q(\theta ,{\theta }^{(i)})=lnL(\theta )=\sum _{k=1}^K(n_{k}ln{a}_k+\sum _{j=1}^N{r}_{jk}ln[\varnothing (y_j|{\theta }_k)]) \\ 约束条件：\sum _{k=1}^K{a}_k=1 \\ 其中n_k=\sum _{j=1}^N{\gamma }_{jk},表示属于第k个分布的样本数 \end{gathered}$$

1.2 高斯混合模型的M步

通过EM算法的E步，我们得到了目标函数$Q(\theta ,{\theta }^{(i)})$，我们要最大化目标函数
$$\begin{gathered} {\theta }^{(i+1)}=argmaxQ(\theta ,{\theta }^{(i)}) \\ s.t.\sum _{k=1}^K{a}_k=1 \\ 求解有约束问题的最优化问题，首先利用拉格朗日乘子法转换为无约束条件 \\ {Q}_0=argmax[Q(\theta ,{\theta }^{(i)})+\lambda (\sum _{k=1}^K{a}_k-1)] \\ {Q}_0=\sum _{k=1}^K(n_{k}ln{a}_k+\sum _{j=1}^N{r}_{jk}ln[\varnothing (y_j|{\theta }_k)])+\lambda (\sum _{k=1}^K{a}_k-1) \\ 将高斯密度函数代入，得到 \\ =\sum _{k=1}^K(n_{k}ln{a}_k+\sum _{j=1}^N{r}_{jk}(-\frac{1}{2}ln2\pi -\frac{1}{2}ln{\sigma }_k^2-\frac{1}{2{\sigma }_k^2}(y_j-u_k^2)))+\lambda (\sum _{k=1}^K{a}_k-1) \end{gathered}$$
然后，我们分别对几个参数求偏导，并令其为0。
$$\begin{gathered} \frac{\partial Q_0}{\partial a_k}=\frac{n_k}{a_k}+\lambda =0 \\ 怎么求\lambda 呢?我们可以利用约束条件\sum _{k=1}^K{a}_k=1。 \\ 易知n_k=-\lambda a_k， \\ 那么\sum _{k=1}^K{n}_k=-\lambda \sum _{k=1}^K{a}_k \\ 左边式子就是总样本数N，右边就变成了-\lambda ，因此-\lambda =N \\ 因此a_k=\frac{n_k}{-\lambda }=\frac{n_k}{N}，和1.2中注下我们计算的一样。 \end{gathered}$$
同理，我们对$u_k$求偏导：
$$\begin{gathered} \frac{\partial Q_0}{\partial u_k}=\sum _{j=1}^N{r}_{jk}(+\frac{2}{2{\sigma }_k^2}(y_j-u_k))=0 \\ 方差{\sigma }^2>0,因此\sum _{j=1}^N{r}_{jk}(y_j-u_k)=0 \\ 因此u_k=\frac{{\sum }_{j=1}^N{r}_{jk}{y}_j}{{\sum }_{j=1}^N{r}_{jk}}=\frac{{\sum }_{j=1}^N{r}_{jk}{y}_j}{n_k}，和1.2中注下我们计算的一样。 \end{gathered}$$
同理，我们对${\sigma }_k^2$求偏导：
$$\begin{gathered} \frac{\partial Q_0}{\partial {\sigma }_k^2}=\sum _{j=1}^N{r}_{jk}(-\frac{1}{2{\sigma }_K^2}+\frac{1}{2{\sigma }_k^4}(y_j-u_k{)}^2)=0 \\ 同理，方差{\sigma }_k^2>0,\sum _{j=1}^N{r}_{jk}(-{\sigma }_k^2+(y_j-u_k{)}^2)=0 \\ 因此{\sigma }_k^2=\frac{{\sum }_{j=1}^N{r}_{jk}(y_j-u_k{)}^2}{{\sum }_{j=1}^N{r}_{jk}}=\frac{{\sum }_{j=1}^N{r}_{jk}(y_j-u_k{)}^2}{n_k}，和1.2中注下我们计算的一样。 \end{gathered}$$

1.3 高斯混合模型举例

我们根据上面公式的推导，现在应该很容易理解下面EM算法求解高斯混合模型参数的过程了。

我们通过一个例子，再重复以下上述求解过程：

假设有6个未知分布的数据 $x1,x2,...,x6$ ，采用2个高斯分布线性组合对其进行拟合。

2 高斯混合模型的应用

1、数据分类： 由于观测数据的分布可以由多个单高斯分布进行近似。通过判断观测数据属于哪个单高斯分布，从而对数据进行分类。

2、异常检测： 数据的异常检测在金融风控、设备状态监测等场景应用广泛。GMM首先对观测数据进行多高斯分布建模，并通过EM算法估计单高斯分布参数，如果测试数据（在线数据）在隶属的高斯分布中的概率密度值小于阈值，则判定为异常点。

3、图像分割：在医学图像分割中应用广泛，通常用于前景以及不同类别的背景之间的分割，本质上与数据分类相同。

4、图像生成： 近期在各大平台上火热的AI绘画，为生成式模型作品。GMM能够用于图像生成关键在于其能够通过不同的高斯分布的线性组合对任意图像分布进行拟合。因此，我们通过GMM求出了图像分布，就可以从图像分布中任意采样生成与原数据相似但不同的图像。

2.1 异常检测

• 我们从正态分布中生成了800个点，从均匀分布中生成了20个点作为异常点。然后为了对比模型效果我们记录了数据是不是异常值的标签(真实数据集一般没有)。

• 准确率99.87%，只有1个点预测错误。而且这个点，和正常值太接近了，预测错了完全可以理解。

• 高斯混合模型是基于分布的，而我们则模拟数据集就是从不同分布里面生成的，所以正好对上了它的专长，所以效果会比较好。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture
from sklearn.metrics import accuracy_score


def get_data():
    # 生成数据：正常数据和异常数据
    np.random.seed(77)
    X_normal = 0.3 * np.random.randn(800, 2)                     # 正常数据
    X_abnormal = np.random.uniform(low=-4, high=4, size=(20, 2)) # 异常数据
    X = np.vstack([X_normal, X_abnormal])

    # 标准化数据
    scaler = StandardScaler()
    scaled = scaler.fit(X_normal)
    X_scaled = scaled.transform(X)
    # 正常点label为1，异常点label为-1
    true_labels = np.hstack([np.ones(len(X_normal)), -1 * np.ones(len(X_abnormal))])
    return X_normal, X_abnormal, X, X_scaled, true_labels



def plot_result(y_pred, model_name=''):
    # Splitting the data into predicted normal and abnormal by gmm
    X_k_normal = X[y_pred == 1]
    X_k_abnormal = X[y_pred == -1]

    # Correctly and incorrectly predicted points by gmm
    correctly_predicted_k = y_pred == true_labels
    incorrectly_predicted_k = ~correctly_predicted_k
    accuracy_k_simple = accuracy_score(true_labels, y_pred)
    print(accuracy_k_simple)

    plt.figure(figsize=(15, 5))
    # Original data plot
    plt.subplot(1, 3, 1)
    plt.scatter(X_normal[:, 0], X_normal[:, 1], color='blue', label='Normal', s=20, marker='x')
    plt.scatter(X_abnormal[:, 0], X_abnormal[:, 1], color='red', label='Abnormal', s=20, marker='x')
    plt.title("Original Data")
    plt.legend()

    # Predicted data plot
    plt.subplot(1, 3, 2)
    plt.scatter(X_k_normal[:, 0], X_k_normal[:, 1], color='b', label='Predicted Normal', s=20, marker='x')
    plt.scatter(X_k_abnormal[:, 0], X_k_abnormal[:, 1], color='r', label='Predicted Abnormal', s=20, marker='x')
    plt.title(f"Predicted Data ({model_name})")
    plt.legend()

    # Prediction accuracy plot
    plt.subplot(1, 3, 3)
    plt.scatter(X[correctly_predicted_k, 0], X[correctly_predicted_k, 1], color='gold', label='Correctly Predicted',
                s=20, marker='x')
    plt.scatter(X[incorrectly_predicted_k, 0], X[incorrectly_predicted_k, 1], color='purple',
                label='Incorrectly Predicted', s=20, marker='x')
    plt.title(f"Prediction Accuracy ({model_name})")
    plt.legend()

    plt.show()


if __name__ == '__main__':
    X_normal, X_abnormal, X, X_scaled, true_labels = get_data()
    # 构建gmm模型进行异常点检测
    gmm = GaussianMixture(n_components=2, random_state=77)
    gmm_labels = gmm.fit_predict(X_scaled)
    plot_result(np.where(gmm_labels == 0, 1, -1), model_name='GaussianMixture')

2.2 数据聚类

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 生成随机数据
np.random.seed(0)
n_samples = 500
X1 = np.random.randn(n_samples, 2)                                # 均值0，方差为1
X2 = np.random.normal(5, 1, 2 * n_samples).reshape(n_samples, 2)  # 均值5，方差为1

X = np.vstack([X1, X2])

# 创建GMM模型对象并拟合数据
n_clusters = 2
gmm = GaussianMixture(n_components=n_clusters).fit(X)

# 预测每个样本所属的聚类
labels = gmm.predict(X)

# 可视化聚类结果
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')
plt.show()

2.3 手写数字生成

• 我们以手写体数字生成为例，利用GMM拟合图像数据分布，并从拟合后的分布中对数据进行采样，从而生成新的手写体数字。
• 我们先下载数据集，并保存本地，方便加载数据。

from sklearn.datasets import fetch_openml
import pickle

# 获取数据集
# False表示以原始格式返回，每个特征是一个单独的数组。True表示返回Pandas
# DataFrame对象
# 自0.24.0（2020 年 12 月）以来，as_frame参数为auto（而不是之前的False默认选项）
mnist = fetch_openml('mnist_784', version=1, as_frame=False)

# 保存数据集到本地
with open('mnist_data.pkl', 'wb') as f:
    pickle.dump(mnist, f)

首先，加载手写体数字数据集，并对图像进行可视化。从70000张数据中选择2000张作为观测数据。

import pickle
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA


def get_data():
    # 从本地导入数据集
    with open('mnist_data.pkl', 'rb') as f:
        mnist = pickle.load(f)
        # 一共7w张图，每个有784个特征。每个特征是28 x 28 像素中的一个点的数值，在0（白）~ 255（黑）之间。
        X, y = mnist["data"], mnist["target"]
    return X, y

def show_images(train_images_part, train_labels_part, picture_type):
    # 显示前10张图片
    fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(8, 4))
    for ax, image, label in zip(axes.ravel(), train_images_part[:10], train_labels_part[:10]):
        ax.imshow(image, cmap=plt.cm.gray_r)
        ax.set_xticks([])
        ax.set_yticks([])
        if label:
            ax.set_title('Label: {}'.format(label))
    plt.suptitle(picture_type, fontsize=16)
    plt.tight_layout()
    plt.show()

if __name__ == '__main__':
    train_images, train_labels = get_data()
    # 1、选取2000张图片
    train_images_part = train_images[:2000]
    train_labels_part = train_labels[:2000]

    # 展示原始图片
    show_images(train_images_part.reshape(-1, 28, 28), train_labels_part, 'original picture')

其次，由于图像为28*28=784维，GMM对高维数据拟合难以收敛，因此，对数据通过PCA进行降维。我们保留85%的主成分，将图像维度降为56维。降维后的图像如下图所示，图像虽部分边缘模糊，但保留了原始图像的主要特征。

    # 2、对图片进行pca降维，将图像维度降为56维
    pca = PCA(0.85, whiten=True)
    train_images_pca = pca.fit_transform(train_images_part)
    # 将经过降维处理的图像转换回原始的高维空间（784维）
    train_images_pca_rec = pca.inverse_transform(train_images_pca)

    show_images(train_images_pca_rec.reshape(-1, 28, 28), train_labels_part, 'pca picture')

最后，我们通过GMM对（2000,56）数据进行拟合，并从中采样10组数据，重建出生成新图像。如下图所示，我们从GMM拟合数据分布中随机采样了10组点，并通过升维重建出了原始图像。可以发现，新生成的图像与原始图像具有较高相似度。

# 3、构建高斯混合模型，并采样10组数据生成新的图片
    gmm = GaussianMixture(110, covariance_type='full', random_state=42)
    gmm.fit(train_images_pca)
    # 采样10组数据生成新的图片
    # data_new[0]为生成的图像
    data_new = gmm.sample(10)
    # 将56维图像转换回原始的高维空间（784维）
    digits_new = pca.inverse_transform(data_new[0])

    show_images(digits_new.reshape(-1, 28, 28), [None] * 10, 'generate picture')