【混合模型】

概述

混合模型（Mixture Model）是一种概率模型，其用于表示总体中子总体的存在，而不需要观测数据识别出该观测数据属于哪一个子总体（子分布）

简单来说，其是一个可以表示观测数据在总体中的概率分布，该概率分布是由 $K$ 个子分布组成的混合分布，对应混合分布的混合模型，就代表了这个总体的概率分布

当需要从子总体的性质推导总体的一些性质时，混合模型能够直接根据总体池的观测值来对子总体的特性进行统计推断，而不需要知道它们的归属信息（属于哪一个子总体）

分层模型

一个典型的有限维度的混合模型是一个分层的模型，其组成如下：

$N$ 个被观测的随机变量，每个随机变量按 $K$ 个子分布构成的混合模型来分布，这些子分布都属于同一类分布，但具体的参数值不同
$N$ 个隐变量，每一隐变量都按 $K$ 维分类分布（隐变量的取值只有 $K$ 个），其说明了对应随机变量所属的子分布是哪一个
$K$ 个参数组，每一参数组都对应一个子分布
$K$ 个混合权重，每个混合权重指定了某个子分布所占的总体的权重，$K$ 个混合权重的和为 $1$

此外，在贝叶斯假设下，混合权重和参数组本身就是随机变量，每个都有一个先验分布，在这种情况下，混合权重可以被视为一个 $K$ 维的随机向量，由分类分布的共轭先验分布得出，参数组将根据各自的先验共轭分布而分布

非贝叶斯假设下的数学表示

从数学角度出发，在非贝叶斯假设下，一个基础的参数化的混合模型可被以下参数所描述：

$K$：混合分布中子分布的个数
$N$：被观测的随机变量的个数
$y_{i=1,\cdots,N}$：第 $i$ 个观测值的随机变量
$z_{i=1,\cdots,N}$：第 $i$ 个观测值所属的子分布
$\boldsymbol{\theta}_{k=1,\cdots,K}$：第 $k$ 个子分布的参数组
$\phi_{k=1,\cdots,K}$：第 $k$ 个子分布的先验概率，$\sum\limits_{k=1}^K \phi_k=1$
$\Phi=[\phi_{1},\phi_{2},\cdots,\phi_{K}]$：由 $\phi_{k}$ 组成的 $K$ 维向量
$F(y|\boldsymbol{\theta})$：某个被观测的随机变量 $y$ 在参数组 $\boldsymbol{\theta}$ 下的概率分布

其中，$y_i|z_i$ 服从 $F(\boldsymbol{\theta}_{z_i})$，即 $y_{i}$ 服从其对应子分布 $z_i$ 的参数组 $\boldsymbol{\theta}_{z_i}$ 指定的概率分布，$z_{i}$ 服从以 $\Phi$ 为概率的 $K$ 类分类分布，即：

$z_{i} \sim \text{Categorical}(\Phi)$

贝叶斯假设下的数学表示

在贝叶斯假设下，所有参数都与随机变量相关，此时，混合模型可被以下参数所描述：

$K$：混合分布中子分布的个数
$N$：被观测的随机变量的个数
$y_{i=1,\cdots,N}$：第 $i$ 个观测值的随机变量
$z_{i=1,\cdots,N}$：第 $i$ 个观测值所属的子分布
$\boldsymbol{\theta}_{k=1,\cdots,K}$：第 $k$ 个子分布的参数组
$\phi_{k=1,\cdots,K}$：第 $k$ 个子分布的先验概率，$\sum\limits_{k=1}^K \phi_k=1$
$\Phi=[\phi_{1},\phi_{2},\cdots,\phi_{K}]$：由 $\phi_{k}$ 组成的 $K$ 维向量
$F(y|\boldsymbol{\theta})$：某个被观测的随机变量 $y$ 在参数组 $\boldsymbol{\theta}$ 下的概率分布
$\alpha$：各子分布参数共用的超参数
$\beta$：混合权重共用的超参数
$H(\boldsymbol{\theta}|\alpha)$：子分布参数组 $\boldsymbol{\theta}$ 的基于 $\alpha$ 的先验概率分布

其中，$\boldsymbol{\theta}$ 服从 $H(\boldsymbol{\theta}|\alpha)$，$y_i|z_i$ 服从 $F(\boldsymbol{\theta}_{z_i})$，即 $y_{i}$ 服从其对应子分布 $z_i$ 的参数组 $\boldsymbol{\theta}_{z_i}$ 指定的概率分布，$z_{i}$ 服从以 $\Phi$ 为概率的 $K$ 类分类分布，即：

$z_{i} \sim \text{Categorical}(\Phi)$

通过使用分布 $F$ 和 $F$ 的共轭先验 $H$，可对观测值和参数进行任意描述，即：

$P(y|\boldsymbol{\theta}) = \sum_{k=1}^K \phi_k F(y|\boldsymbol{\theta}_k)$

对离散观测值和连续观测值来说，常见的 $F$ 是：

离散观测值：分类分布
连续观测值：高斯分布

此外，二项分布、多项分布、泊松分布、指数分布等均可作为 $F$

【高斯混合模型】

概述

高斯混合模型（Gaussian mixture model，GMM）是混合模型的一种，可用于拟合未知的参数向量或多元正态分布，目前常用于聚类，其通过概率模型来表达聚类原型，因此去其是原型聚类的一种

其与 K-Means 的相同之处在于，都需要指定 K 值，但其通过 EM 算法进行参数估计，有时只能收敛到局部最优

此外，与 K-Means 不同，GMM 是通过选择成分最大化后验概率来完成聚类的，各数据点的后验概率表示属于各类的可能性，而不是判定它完全属于某个类，所以其是软聚类

在各类尺寸不同、聚类间有相关关系的时候，GMM 要比 K-Means 更为合适

概率分布

GMM 是一个生成模型，其假设数据是从 $K$ 个高斯分布中生成的，为每一个高斯分布赋予一个权重，每当生成一个数据时，就按权重的比例随机选择一个分布，然后按照该分布生成数据

根据数据进行反推，假设一个样本有一个潜在的类别，这个类别是无法观测到的隐变量，那么对于给定样本的观测值 $y$，在给定参数组 $\boldsymbol{\theta}=(\mu,\sigma^2)$ 的条件下，边际概率为：

$\begin{align*} P(y|\Theta) &= \sum_{k}^K P(y,z=C_k|\boldsymbol{\theta}_k) \\ &= \sum_{k}^K P(z=C_k|\boldsymbol{\theta}_k) P(y|\boldsymbol{\theta}_k,z=C_k) \end{align*}$

其中，$z=C_k$ 表示样本所属的类别，$P(z=C_k|\boldsymbol{\theta}_k)$ 即隐变量的先验概率 $\phi_k$

由于样本的条件概率分布服从于高斯概率分布，那么对于 $P(y|\boldsymbol{\theta}_k,z=C_k)$，有：

$P(y|\boldsymbol{\theta}_k,z=C_k) = \frac{1}{\sqrt{2\pi} \sigma_k} \exp \Big[ -\frac{(y-\mu_k)^2}{2\sigma_k^2} \Big]$

将其记为 $f(y|\boldsymbol{\theta}_k)$，则样本的边际概率为：

$P(y|\boldsymbol{\Theta}) = \sum_{k=1}^K \phi_k f(y|\boldsymbol{\theta}_k)$

其中，$\Theta=(\phi_1,\phi_2,\cdots,\phi_k;\boldsymbol{\theta}_1,\boldsymbol{\theta}_2,\cdots,\boldsymbol{\theta}_k)$，$\phi_k$ 是第 $k$ 个子分布的先验概率，$f(y|\boldsymbol{\theta}_k)$ 是被观测的随机变量 $y$ 在参数组 $\boldsymbol{\theta}_k$ 下的高斯概率密度函数，其由第 $k$ 个子分布的参数组 $\boldsymbol{\theta}_k=(\mu_k,\sigma^2_k)$ 来决定

【高斯混合模型的参数估计】

概述

假设观测数据 $y_1,y_2,\cdots,y_n$ 由高斯混合模型生成

$P(y|\boldsymbol{\Theta}) = \sum_{k=1}^K \phi_k f(y|\boldsymbol{\theta}_k)$

对于单高斯分布模型，可以使用最大似然法来估算参数 $\theta$ 的值，但对于高斯混合模型来说，由于事先并不知道每个观测点属于哪个子分布，同时每个子分布中都有未知的 $\phi_k,\mu_k,\sigma_k^2$，因此无法像单高斯分布模型那样使用最大似然法来求导来估算参数 $\Theta$，需要通过迭代的方式求解

目前，对高斯混合模型参数估计常采用 EM 算法

基本思想

确定对数似然函数

对于观测数据 $y_i,i=1,2,\cdots,n$，假设其产生流程为：依概率 $\phi_k$ 选择第 $k$ 个高斯分布模型 $f(y|\theta_k)$，然后依第 $k$ 个分模型的概率分布 $f(y|\theta_k)$ 生成观测数据 $y_i$

此时，观测数据 $y_i$ 是已知的，反映观测数据 $y_i$ 来自第 $k$ 个分模型的数据是未知的，以隐变量 $z_{ik}$ 来表示，有：

$z_{ik} = \left\{\begin{array}{cc} 1,& 第 i 个观测来自第k个分模型\\ 0,& 否则 \end{array}\right.$

其中，隐变量 $z_{ik},i=1,2,\cdots,n,k=1,2,\cdots,K$ 是 0-1 随机变量

那么此时，对于观测数据 $y_i$ 和隐变量 $z_{ik}$，完全数据为：

$(y_i,z_{i1},z_{i2},\cdots,z_{iK}),\quad i=1,2,\cdots,n$

于是，可以写出完全数据的似然函数：

$\begin{align*} L(y,z|\Theta) &= \prod_{i=1}^n P(y_i,z_{i1},z_{i2},\cdots,z_{iK}|\Theta) \\ &= \prod_{k=1}^K \prod_{i=1}^n [\phi_k f(y_i|\boldsymbol{\theta}_{k})]^{z_{ik}} \\ &= \prod_{k=1}^K \phi_{k}^{n_k} \prod_{i=1}^n [f(y_i|\boldsymbol{\theta}_{k})]^{z_{ik}} \\ &= \prod_{k=1}^K \phi_{k}^{n_k} \prod_{i=1}^n \Bigg[ \frac{1}{\sqrt{2\pi} \sigma_k} \exp \Big[ -\frac{(y-\mu_k)^2}{2\sigma_k^2} \Big] \Bigg]^{z_{ik}} \end{align*}$

其中，对于 $n_k$ ，有：

$\begin{gather*} n_k = \sum\limits_{i=1}^n z_{ik} \\ \sum_{k=1}^K n_k = n \end{gather*}$

此时完全数据的对数似然函数为：

$\ln L(y,z|\Theta) = \sum_{k=1}^K \Bigg\{ n_k \ln \phi_k + \sum_{i=1}^n z_{ik} \Big[ \ln \frac{1}{\sqrt{2\pi}}-\ln \sigma_k -\frac{(y_i-\mu_k)^2}{2\sigma_k^2} \Big]\Bigg\}$

E 步：确定 Q 函数

在得出完全数据的对数似然函数后，可得到 Q 函数，即：

$\begin{align*} Q(\Theta,\Theta_j) &= E[\ln L(y,z|\Theta) | y,\Theta_j] \\ &= E\Bigg\{ \sum_{k=1}^K \Big\{ n_k \ln \phi_k + \sum_{i=1}^n z_{ik}\big[ \ln\frac{1}{\sqrt{2\pi}} - \ln \sigma_k -\frac{(y_i-\mu_k)^2}{2\sigma_k^2} \big] \Big\} \Bigg\} \\ &= \sum_{k=1}^K \Big\{ \sum_{i=1}^n Ez_{ik} \ln \phi_k + \sum_{i=1}^n Ez_{ik} \big[ \ln\frac{1}{\sqrt{2\pi}} - \ln \sigma_k -\frac{(y_i-\mu_k)^2}{2\sigma_k^2} \big]\Big\} \\ \end{align*}$

此时，对于隐变量 $z_{ik}$ 的期望 $Ez_{ik}$，有：

$\begin{align*} Ez_{ik} &= E(z_{ik}|y,\Theta) \\ &= P(z_{ik}=1|y,\Theta) \\ &= \frac{P(z_{ik}=1,y_i|\Theta)}{\sum_{k=1}^K \limits P(z_{ik}=1,y_i|\Theta)} \\ &= \frac{P(y_i|z_{ik}=1,\Theta)P(z_{ik}=1|\Theta)}{\sum_{k=1}^K \limits P(y_i|z_{ik}=1,\Theta)P(z_{ik}=1|\Theta)} \\ &= \frac{\phi_k f(y_i|\theta_k)}{\sum\limits_{k=1}^K \phi_k f(y_i|\theta_k)} \end{align*}$

将其记为 $\hat{z}_{ik}$，其是在当前模型参数下第 $i$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $y_i$ 的响应度

将 $\hat{z}_{ik}=Ez_{ik}$ 与 $n_k =\sum\limits_{i=1}^n Ez_{ik}$ 代回 $Q(\Theta,\Theta_j)$，即有：

$Q(\Theta,\Theta_j) = \sum_{k=1}^K \Bigg\{ n_k \ln \phi_k + \sum_{i=1}^n \hat{z}_{ik} \Big[ \ln \frac{1}{\sqrt{2\pi}} - \ln \sigma_k -\frac{(y_i-\mu_k)^2}{2\sigma_k^2} \Big] \Bigg\}$

M 步：求极大值

迭代的 M 步是求 Q 函数 $Q(\Theta,\Theta_j)$ 对 $\Theta$ 的极大值，即对每一轮迭代的模型参数

$\Theta_{j+1} = \arg \max_{\Theta} Q(\Theta,\Theta_j)$

对于 Q 函数，分别用 $\hat{\phi}_k,\hat{\mu}_{k},\hat{\sigma}_k^2,k=1,2,\cdots,K$ 来表示 $\Theta_{j+1}$ 的各参数，令它们可以通过对 Q 函数求偏导，并令结果为 $0$ 来得到

具体来说，在 $\sum\limits_{k=1}^K \phi_k = 1$ 条件下，求偏导数并令其等于 $0$，有：

$\hat{\phi}_{k} = \frac{n_k}{n} = \frac{1}{n} \sum_{i=1}^n \hat{z}_{ik}$

分别对 $\mu_k,\sigma_k^2$ 求偏导数并令其等于 $0$，有：

$\begin{gather*} \hat{\mu}_k = \frac{\sum\limits_{i=1}^n \hat{z_{ik}}y_i}{\sum\limits_{i=1}^n \hat{z}_{ik}} \\ \hat{\sigma}_k^2 = \frac{\sum\limits_{i=1}^n \hat{z_{ik}} (y_i-\mu_k)^2}{\sum\limits_{i=1}^n \hat{z}_{ik}} \end{gather*}$

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

算法流程

输入：观测数据 $y_1,y_2,\cdots,y_n$，高斯混合模型

$\begin{align*} P(y|\boldsymbol{\Theta}) &= \sum\limits_{k=1}^K \phi_k f(y|\boldsymbol{\theta}_k) \\ \boldsymbol{\theta}_k&=(\mu_k,\sigma^2_k),\quad k=1,2,\cdots,K \end{align*}$

输出：高斯混合模型参数 $\Theta=(\phi_1,\phi_2,\cdots,\phi_k;\boldsymbol{\theta}_1,\boldsymbol{\theta}_2,\cdots,\boldsymbol{\theta}_k)$

算法步骤：

Step1：随机初始化模型参数 $\Theta$ 的初值 $\Theta_0$

Step2：E 步，依据当前模型参数，计算分模型 $k$ 对观测数据 $y_i$ 的响应度

$\hat{z}_{ik}=\frac{\phi_k f(y_i|\theta_k)}{\sum\limits_{k=1}^K \phi_k f(y_i|\theta_k)},\quad i=1,2,\cdots,n,k=1,2,\cdots,K$

Step3：M 步，根据响应度 $\hat{z}_{ik}$ 计算新一轮迭代的模型参数 $\hat{\phi}_k,\hat{\mu}_{k},\hat{\sigma}_k^2,k=1,2,\cdots,K$

$\begin{gather*} \hat{\phi}_{k} = \frac{1}{n} \sum_{i=1}^n \hat{z}_{ik} \\ \hat{\mu}_k = \frac{\sum\limits_{i=1}^n \hat{z_{ik}}y_i}{\sum\limits_{i=1}^n \hat{z}_{ik}} \\ \hat{\sigma}_k^2 = \frac{\sum\limits_{i=1}^n \hat{z_{ik}} (y_i-\mu_k)^2}{\sum\limits_{i=1}^n \hat{z}_{ik}} \end{gather*}$

Step4：重复 Step2、Step3，直到收敛

【sklearn 实现】

以 sklearn 中的鸢尾花数据集为例，选取其后两个特征来实现利用 GMM 聚类

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture as GMM
from scipy.spatial.distance import cdist
from matplotlib.colors import ListedColormap

# 特征提取
def deal_data():
    iris = load_iris()  # sklearn的鸢尾花数据集
    # iris分为三类，前50行一类，51-100行一类，101-150行一类
    X = iris.data[:, [2, 3]] # 选用后两个特征作为样本特征
    return X

# 数据标准化
def standard_scaler(X):
    sc = StandardScaler() # 初始化一个sc对象去对数据集作变换
    scaler = sc.fit(X) # 归一化，存有计算出的均值和方差
    X_std = scaler.transform(X) # 利用 scaler 进行标准化
    return X_std

# 模型训练
def train_model(features, k):
    # 建立聚类器
    model = GMM(n_components=k)
    # 聚类
    model.fit(features)
    return model

# 可视化
def visualization(X, y):
    # 创建 color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

     # 全数据集，不同类别样本点的特征作为坐标(x,y)，用不同颜色画散点图
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=cmap(idx), marker=markers[idx], label=cl) 
        
    plt.xlabel('petal length [standardized]')
    plt.ylabel('petal width [standardized]')
    plt.legend(loc='upper left')
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    # 数据处理
    X = deal_data()
    
    # 数据标准化
    X_std = standard_scaler(X)

    # 模型训练
    k = 3
    model = train_model(X_std, k)

    # 获取聚类标签
    label_pred = model.predict(X_std) 
    
    # 打印前10个点分属于每个类的概率
    probs = model.predict_proba(X_std)
    print(probs[:10].round(2))
    
    # 可视化
    visualization(X_std, label_pred)

Alex_McAvoy

高斯混合模型 GMM