24/11/14 算法笔记 GMM高斯混合模型

高斯混合模型（Gaussian Mixture Model，简称 GMM）是一种概率模型，用于表示具有多个子群体的数据集，其中每个子群体的数据分布可以用高斯分布（正态分布）来描述。GMM 是一种软聚类方法，意味着它为每个数据点分配一个属于每个聚类的概率分布，而不是硬聚类方法中的严格分类。

GMM 的组成

一个 GMM 由以下几个部分组成：

1. 聚类数量（K）：模型中高斯分布（聚类）的数量。
2. 均值向量（μkμk​）：每个高斯分布的均值向量，其中 kk 表示聚类索引。
3. 协方差矩阵（ΣkΣk​）：每个高斯分布的协方差矩阵，描述了数据在各个维度上的分布范围和形状。
4. 混合系数（πkπk​）：每个高斯分布的权重，表示数据属于该聚类的概率，所有混合系数之和为1。

GMM 的数学表达

GMM 的概率密度函数（PDF）可以表示为：

GMM 的学习

GMM 的参数学习通常使用 EM 算法进行，EM算法前面有将，是一个策略优化算法

24/11/14 算法笔记 EM算法期望最大化算法-CSDN博客

我们来看一下简单的GMM源代码

import numpy as np
from scipy.stats import multivariate_normal

class GaussianMixture:
    def __init__(self, n_components, covariance_type='full', n_iter=100, random_state=None):
        self.n_components = n_components  # 聚类数量
        self.covariance_type = covariance_type  # 协方差类型
        self.n_iter = n_iter  # 迭代次数
        self.random_state = random_state  # 随机种子
        self.weights_ = None  # 混合系数
        self.means_ = None  # 均值
        self.covariances_ = None  # 协方差

    def _initialize_parameters(self, X):
        """随机初始化均值、协方差和权重"""
        n_samples, n_features = X.shape
        self.weights_ = np.ones(self.n_components) / self.n_components  # 初始化权重
        random_indices = np.random.choice(n_samples, self.n_components, replace=False)
        self.means_ = X[random_indices]  # 随机选择均值
        self.covariances_ = np.array([np.eye(n_features)] * self.n_components)  # 初始化协方差为单位矩阵

    def _e_step(self, X):
        """E步骤：计算每个数据点属于每个高斯分布的责任"""
        n_samples = X.shape[0]
        responsibilities = np.zeros((n_samples, self.n_components))
        
        for k in range(self.n_components):
            rv = multivariate_normal(mean=self.means_[k], cov=self.covariances_[k])
            responsibilities[:, k] = self.weights_[k] * rv.pdf(X)
        
        # 归一化责任
        responsibilities /= responsibilities.sum(axis=1, keepdims=True)
        return responsibilities

    def _m_step(self, X, responsibilities):
        """M步骤：更新均值、协方差和权重"""
        n_samples = X.shape[0]
        effective_n = responsibilities.sum(axis=0)  # 每个聚类的有效样本数量
        
        # 更新权重
        self.weights_ = effective_n / n_samples
        
        # 更新均值
        self.means_ = np.dot(responsibilities.T, X) / effective_n[:, np.newaxis]
        
        # 更新协方差
        for k in range(self.n_components):
            diff = X - self.means_[k]
            self.covariances_[k] = np.dot(responsibilities[:, k] * diff.T, diff) / effective_n[k]

    def fit(self, X):
        """训练模型"""
        self._initialize_parameters(X)  # 初始化参数
        for _ in range(self.n_iter):  # 迭代更新
            responsibilities = self._e_step(X)  # E步骤
            self._m_step(X, responsibilities)  # M步骤

    def predict(self, X):
        """预测数据点的聚类标签"""
        responsibilities = self._e_step(X)  # 计算责任
        return np.argmax(responsibilities, axis=1)  # 返回最大责任的聚类索引

    def sample(self, n_samples):
        """从模型中生成新样本"""
        samples = np.zeros((n_samples, self.means_.shape[1]))
        for i in range(n_samples):
            k = np.random.choice(self.n_components, p=self.weights_)  # 根据权重选择聚类
            samples[i] = np.random.multivariate_normal(self.means_[k], self.covariances_[k])  # 生成样本
        return samples

接下来让我们分析下每段代码

1.初始化函数 __init__

def __init__(self, n_components, covariance_type='full', n_iter=100, random_state=None):
    self.n_components = n_components  # 聚类数量
    self.covariance_type = covariance_type  # 协方差类型
    self.n_iter = n_iter  # 迭代次数
    self.random_state = random_state  # 随机种子
    self.weights_ = None  # 混合系数
    self.means_ = None  # 均值
    self.covariances_ = None  # 协方差

这是类的构造函数，用于初始化GMM模型的参数：

• n_components：模型中高斯分布（聚类）的数量。
• covariance_type：协方差矩阵的类型，可以是'full'、'diag'或'spherical'，分别表示全协方差、对角协方差和球面协方差。
• n_iter：EM算法的最大迭代次数。
• random_state：随机数生成器的种子，用于结果的可重复性。
• weights_、means_和covariances_：这些属性将在模型训练后存储模型参数。

2.参数初始化函数 _initialize_parameters

def _initialize_parameters(self, X):
    """随机初始化均值、协方差和权重"""
    n_samples, n_features = X.shape
    self.weights_ = np.ones(self.n_components) / self.n_components  # 初始化权重
    random_indices = np.random.choice(n_samples, self.n_components, replace=False)
    self.means_ = X[random_indices]  # 随机选择均值
    self.covariances_ = np.array([np.eye(n_features)] * self.n_components)  # 初始化协方差为单位矩阵

这个函数用于随机初始化模型参数：

• self.weights_：权重初始化为均等分布。
• self.means_：均值初始化为数据集中随机选择的点。
• self.covariances_：协方差矩阵初始化为单位矩阵，适用于全协方差情况。
• 协方差可以告诉我们两个变量是如何一起变化的。如果两个变量的协方差是正的，那么它们倾向于朝相同的方向变化；如果协方差是负的，那么一个变量增加时，另一个变量倾向于减少。

3.E步骤函数 _e_step

def _e_step(self, X):
    """E步骤：计算每个数据点属于每个高斯分布的责任"""
    n_samples = X.shape[0]
    responsibilities = np.zeros((n_samples, self.n_components))
    
    for k in range(self.n_components):
        #函数用于生成符合多元正态分布的随机样本。
        rv = multivariate_normal(mean=self.means_[k], cov=self.covariances_[k])
        responsibilities[:, k] = self.weights_[k] * rv.pdf(X)
    
    # 归一化责任
    responsibilities /= responsibilities.sum(axis=1, keepdims=True)
    return responsibilities

E步骤计算每个数据点属于每个高斯分布的责任（后验概率）：

• 使用multivariate_normal.pdf计算每个高斯分布的PDF值。
• 将每个高斯分布的PDF值乘以相应的权重，得到未归一化的责任。
• 通过将每个数据点的责任除以其总和来归一化责任，确保每个数据点的责任之和为1。

PDF值通常指的是概率密度函数（Probability Density Function）的值。概率密度函数是连续概率分布的一个核心概念，它描述了随机变量在给定区间内取值的概率密度。对于连续随机变量，其概率密度函数的图形可以告诉我们随机变量取某个特定值的可能性。

4.M步骤函数 _m_step

def _m_step(self, X, responsibilities):
    """M步骤：更新均值、协方差和权重"""
    n_samples = X.shape[0]
    effective_n = responsibilities.sum(axis=0)  # 每个聚类的有效样本数量
    
    # 更新权重
    self.weights_ = effective_n / n_samples
    
    # 更新均值
    self.means_ = np.dot(responsibilities.T, X) / effective_n[:, np.newaxis]
    
    # 更新协方差
    for k in range(self.n_components):
        diff = X - self.means_[k]
        self.covariances_[k] = np.dot(responsibilities[:, k] * diff.T, diff) / effective_n[k]

M步骤根据E步骤计算的责任更新模型参数：

• self.weights_：权重更新为每个聚类的有效样本数量除以总样本数量。
• self.means_：均值更新为加权平均，权重是每个数据点对每个聚类的责任。
• self.covariances_：协方差更新为加权的样本偏差的外积，权重是每个数据点对每个聚类的责任。

5.训练函数 fit

def fit(self, X):
    """训练模型"""
    self._initialize_parameters(X)  # 初始化参数
    for _ in range(self.n_iter):  # 迭代更新
        responsibilities = self._e_step(X)  # E步骤
        self._m_step(X, responsibilities)  # M步骤

• 首先调用_initialize_parameters函数初始化参数。
• 然后进行指定次数的迭代，每次迭代都包括E步骤和M步骤。

6.预测函数 predict

def predict(self, X):
    """预测数据点的聚类标签"""
    responsibilities = self._e_step(X)  # 计算责任
    return np.argmax(responsibilities, axis=1)  # 返回最大责任的聚类索引

• 首先调用_e_step函数计算新数据点对每个聚类的责任。
• 然后返回责任最大的聚类索引作为预测标签。

7.采样函数 sample

def sample(self, n_samples):
    """从模型中生成新样本"""
    samples = np.zeros((n_samples, self.means_.shape[1]))
    for i in range(n_samples):
        k = np.random.choice(self.n_components, p=self.weights_)  # 根据权重选择聚类
        samples[i] = np.random.multivariate_normal(self.means_[k], self.covariances_[k])  # 生成样本
    return samples

• 首先初始化一个空的样本数组。
• 然后根据每个聚类的权重随机选择一个聚类。
• 从选定的聚类对应的高斯分布中生成一个样本。
• 重复上述过程，直到生成所需数量的样本。