高斯混合模型 GMM计算方法

高斯混合模型（Gaussian Mixture Model, GMM）是一种常用的概率模型，用于聚类和密度估计。它假设数据是由多个高斯分布混合生成的。GMM 的计算通常使用 期望最大化（Expectation-Maximization, EM）算法 来求解。

一、问题设定

我们有：

数据集：\(X = \{x_1, x_2, ..., x_N\}\)，其中每个$ x_i \in \mathbb{R}^D$

假设有 \(K\) 个高斯分布

每个样本由一个隐变量 \(z_i \in \{1,2,...,K\}\) 表示其来自哪个高斯分布
我们的目标是：

• 推断出所有高斯分布的参数：均值 \(\mu_k\) 、协方差矩阵 \(\Sigma_k\) 、权重 \(\alpha_k\)

• 同时推断每个样本属于哪一个高斯分布的概率

二、GMM 的联合分布形式

GMM 的联合分布为：

\[p(x) = \sum_{k=1}^K \alpha_k \cdot \mathcal{N}(x | \mu_k, \Sigma_k) \]

其中：

• \(\alpha_k\) 是第\(k\) 个高斯分布的权重，满足 \(\sum_{k=1}^K \alpha_k = 1\)

• \(\mathcal{N}(x | \mu_k, \Sigma_k)\) 是多元高斯分布：

\[\mathcal{N}(x | \mu_k, \Sigma_k) = \frac{1}{(2\pi)^{D/2} |\Sigma_k|^{1/2}} \exp\left(-\frac{1}{2}(x - \mu_k)^T \Sigma_k^{-1} (x - \mu_k)\right) \]

三、EM 算法流程（E 步 + M 步）

EM 算法是一个迭代优化算法，用于在存在隐变量的情况下进行最大似然估计。

✅ 初始化参数（Iteration 0）
随机初始化以下参数：

• 高斯权重：\(\alpha_k^{(0)}\) ，满足 \(\sum \alpha_k = 1\)
• 均值向量：\(\mu_k^{(0)} \in \mathbb{R}^D\)
• 协方差矩阵：\(\Sigma_k^{(0)} \in \mathbb{R}^{D \times D}\) ，正定对称矩阵（可初始化为单位矩阵）

四、E 步：计算后验概率（责任函数）

对于每个样本 \(x_i\) ，计算其属于第$ k$ 个高斯分布的后验概率：

\[\gamma(z_{ik}) = p(z_i = k | x_i) = \frac{\alpha_k \cdot \mathcal{N}(x_i | \mu_k, \Sigma_k)}{\sum_{j=1}^K \alpha_j \cdot \mathcal{N}(x_i | \mu_j, \Sigma_j)} \]

这个值表示：第 \(i\) 个样本“属于”第 \(k\) 个高斯分布的“责任”。

五、M 步：更新参数

定义：

• \(N_k = \sum_{i=1}^N \gamma(z_{ik})\)：第 \(k\) 个高斯分布的“有效样本数”
  更新公式如下：

1. 更新权重：

\[ \alpha_k^{(t+1)} = \frac{N_k}{N} \]

2. 更新均值：

\[\mu_k^{(t+1)} = \frac{1}{N_k} \sum_{i=1}^N \gamma(z_{ik}) x_i \]

3. 更新协方差矩阵：

\[\Sigma_k^{(t+1)} = \frac{1}{N_k} \sum_{i=1}^N \gamma(z_{ik}) (x_i - \mu_k)(x_i - \mu_k)^T \]

六、迭代与终止条件

重复 E 步 和 M 步，直到满足以下任意一种情况：

• 参数变化小于某个阈值（例如：\(||\theta^{(t+1)} - \theta^{(t)}|| < \epsilon\) )
• 对数似然变化很小（如：\(\log p(X|\theta^{(t+1)}) - \log p(X|\theta^{(t)}) < \epsilon\) ）
• 达到预设的最大迭代次数

七、对数似然函数（可用于评估收敛）

每一步可以计算当前参数下的对数似然：

\[\log p(X | \theta) = \sum_{i=1}^N \log \left( \sum_{k=1}^K \alpha_k \cdot \mathcal{N}(x_i | \mu_k, \Sigma_k) \right) \]

这个值应随着迭代逐渐上升并趋于稳定。

八、最终输出

经过多次迭代后，我们得到一组最优参数估计：

• 最终的权重：\(\alpha_1, \alpha_2, ..., \alpha_K\)
• 最终的均值：\(\mu_1, \mu_2, ..., \mu_K\)
• 最终的协方差矩阵：\(\Sigma_1, \Sigma_2, ..., \Sigma_K\)

同时我们可以将每个样本分配给最可能的类别（即选择 $ \arg\max_k \gamma(z_{ik})$。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from matplotlib.patches import Ellipse
# 高斯分布概率密度函数
def gaussian_pdf(X, mu, Sigma):
    return multivariate_normal.pdf(X, mean=mu, cov=Sigma)
# E-step
def e_step(X, alpha, mu, Sigma):
    N, K = X.shape[0], len(alpha)
    gamma = np.zeros((N, K))
    for i in range(N):
        total = 0.0
        for k in range(K):
            gamma[i, k] = alpha[k] * gaussian_pdf(X[i], mu[k], Sigma[k])
            total += gamma[i, k]
        gamma[i, :] /= total
    return gamma

# M-step
def m_step(X, gamma):
    N, D = X.shape
    K = gamma.shape[1]
    alpha_new = np.zeros(K)
    mu_new = np.zeros((K, D))
    Sigma_new = [np.zeros((D, D)) for _ in range(K)]
    for k in range(K):
        Nk = gamma[:, k].sum()
        alpha_new[k] = Nk / N
        mu_new[k] = np.dot(gamma[:, k], X) / Nk
        diff = X - mu_new[k]
        Sigma_new[k] = np.dot(gamma[:, k] * diff.T, diff) / Nk
    return alpha_new, mu_new, Sigma_new

# 绘图函数
def plot_gmm(X, mu, Sigma, labels=None, title=""):
    plt.figure(figsize=(8, 6))
    if labels is not None:
        colors = ['r', 'g', 'b', 'c', 'm', 'y', 'k']
        for k in range(len(mu)):
            idx = (labels == k)
            plt.scatter(X[idx, 0], X[idx, 1], c=colors[k], label=f"Cluster {k+1}", s=50)
    else:
        plt.scatter(X[:, 0], X[:, 1], c='gray', s=30)
    for k in range(len(mu)):
        plot_cov_ellipse(Sigma[k], mu[k], nstd=2, alpha=0.3)
    plt.title(title)
    plt.xlabel("Feature 1")
    plt.ylabel("Feature 2")
    plt.legend()
    plt.grid(True)
    plt.pause(0.5)

def plot_cov_ellipse(cov, pos, nstd=2, color='blue', **kwargs):
    eigvals, eigvecs = np.linalg.eigh(cov)
    order = eigvals.argsort()[::-1]
    eigvals, eigvecs = eigvals[order], eigvecs[:, order]
    theta = np.degrees(np.arctan2(*eigvecs[:, 0][::-1]))
    width, height = 2 * nstd * np.sqrt(eigvals)
    ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, color=color, **kwargs)
    plt.gca().add_patch(ellip)


def main():

    # 数据集
    X = np.array([
        [1.0, 2.0],
        [1.5, 1.8],
        [5.0, 8.0],
        [8.0, 8.0],
        [1.0, 0.6],
        [9.0, 11.0]
    ])

    N, D = X.shape
    K = 2

    # 初始化参数
    np.random.seed(42)
    alpha = np.random.rand(K)
    alpha /= alpha.sum()
    indices = np.random.choice(N, K, replace=False)
    mu = X[indices]
    Sigma = [np.eye(D) for _ in range(K)]

    # 主训练循环 + 可视化
    max_iter = 20
    tolerance = 1e-4
    log_likelihood_prev = -np.inf

    for it in range(max_iter):
        gamma = e_step(X, alpha, mu, Sigma)
        alpha, mu, Sigma = m_step(X, gamma)

        log_likelihood = 0.0
        for i in range(N):
            ll = 0.0
            for k in range(K):
                ll += alpha[k] * gaussian_pdf(X[i], mu[k], Sigma[k])
            log_likelihood += np.log(ll)

        print(f"Iteration {it+1}, Log Likelihood: {log_likelihood:.4f}")
        labels = np.argmax(gamma, axis=1)
        plot_gmm(X, mu, Sigma, labels=labels, title=f"Iteration {it+1}")

        if abs(log_likelihood - log_likelihood_prev) < tolerance:
            print("Converged.")
            break

        log_likelihood_prev = log_likelihood

    plt.show()

if __name__ == '__main__':
    main()