用概率密度建模的高斯混合模型

可视化演示

展示GMM用多个高斯分布混合建模数据，通过EM算法迭代估计各分量参数实现软聚类。

高斯分布基石

舞台上的小提琴组：二十位乐手在同一段落里演奏，虽然每个人的音准都有细微偏差，但大部分音符密集地围绕在同一个主音附近。把这些音符画成柱状图，就会得到一座中间高、两边矮的“钟形山”。这座山的主人——高斯分布——是整部“数据交响曲”里最小的演奏单位，也是 GMM 的起点。

在数学上，描述这座钟形山的工具叫概率密度函数。山上任意位置的高度，代表“音符落在这里的可能性有多密”。山峰的位置叫均值 μ（主音高），山的胖瘦叫方差 σ²（音准偏差有多大）。形状陡峭说明乐手们很齐，矮胖说明大家各吹各的调。

$$p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

1. 概率密度 $p(x)$：输出的是该点附近区域的“密度”，数值越大代表数据越密集地聚集在这里。

$$p(x)$$

2. 均值 $\mu$：山峰的位置，也就是这个声部的“主音高”，所有数据围绕它散布。

$$\mu$$

3. 方差 $\sigma^2$：山的胖瘦，刻画乐手们音准偏差的大小，方差越大钟形越矮胖。

$$\sigma^2$$

4. 归一化系数 $\frac{1}{\sqrt{2\pi\sigma^2}}$：保证整座山包住的面积恒为 1，使 $p(x)$ 成为合法的概率密度。

$$\frac{1}{\sqrt{2\pi\sigma^2}}$$

5. 指数项 $e^{-\frac{(x-\mu)^2}{2\sigma^2}}$：离均值越远，函数值衰减得越快，正是它塑造了钟形轮廓。

$$e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

当数据有多个维度时（比如每个用户的两个属性），方差会升级为协方差矩阵 Σ，钟形山变成“钟形面团”，可以被压扁、拉长、旋转——这正是椭圆簇的来源。一支小提琴组就是一个高斯分布：音高定在哪、乐手们多默契，决定这座山的模样。而一场完整的演出，自然是由许多支这样的小组叠加出来的。

混合模型思想

交响乐从来不是只有小提琴。长笛、圆号、定音鼓各占一个声部，调音台上的三个音量推子，决定每个声部在总响度里占几成。整体听感 = 各声部声音 × 各自音量。GMM 的“混合”二字，就是把这种混音逻辑写成数学公式。

GMM 说：数据集里其实藏着 K 组“乐器”（高斯成分），每个数据点就像从其中某组乐器里“演奏”出来的音符。第 k 组的音量称为混合权重 π_k，权重越大，从它那里冒出来的数据点就越多。把所有声部叠起来，就得到整个数据集的地图。

$$p(x) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(x \mid \mu_k, \Sigma_k), \qquad \sum_{k=1}^{K}\pi_k = 1, \quad \pi_k \ge 0$$

1. 成分数量 $K$：台上乐器组的总数，也就是我们要聚类的簇数。

$$K$$

2. 混合权重 $\pi_k$：第 $k$ 组乐器的音量推子，表示该成分在全局中占的比例。

$$\pi_k$$

3. 高斯成分 $\mathcal{N}(x \mid \mu_k, \Sigma_k)$：第 $k$ 组乐器自己的钟形山，多维时用协方差矩阵 $\Sigma_k$ 控制形状。

$$\mathcal{N}(x \mid \mu_k, \Sigma_k)$$

4. 求和叠加 $\sum_{k=1}^{K}$：把 K 座山按各自权重叠成一座连绵的群山。

$$\sum_{k=1}^{K}$$

5. 约束条件 $\sum_{k=1}^{K}\pi_k = 1, \ \pi_k \ge 0$：所有推子加起来正好是总音量，并且不允许出现“负音量”。

$$\sum_{k=1}^{K}\pi_k = 1, \quad \pi_k \ge 0$$

调音台上的推子就是权重 π：推子拉满的声部贡献大，推子几乎归零的声部近乎无声。GMM 的模型写下来，不过是把“混音逻辑”记录在案。但真正的挑战来了——你手上只有混好的录音，怎么反推出每个推子的位置？这就轮到软聚类与 EM 登场了。

软聚类原理

现在反过来：你只有混好的录音，要判断某个音符是哪件乐器发出的。贝斯手按下一个揉弦长音，低音提琴可能也在同频附和——你很难说“这 100% 是贝斯”，更自然的说法是“70% 贝斯、30% 低音提琴”。

每个数据点 x_i 对每个成分 k 都有一个“归属度”，数学上叫后验概率 γ_ik，常被称作“责任度”——第 k 个成分对生成这个点负多大责任。这就是软聚类：一个点可以同时属于多个簇，只是概率有大小之分。硬聚类只允许选一个乐器，软聚类则保留“兼收并蓄”的信息。

$$\gamma_{ik} = \frac{\pi_k \, \mathcal{N}(x_i \mid \mu_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(x_i \mid \mu_j, \Sigma_j)}$$

1. 责任度 $\gamma_{ik}$：第 $i$ 个点对第 $k$ 个簇的归属概率，值在 0 到 1 之间。

$$\gamma_{ik}$$

2. 分子 $\pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$：这个声部的“吸引力”——音量越大、音符离它越近，贡献越大。

$$\pi_k \, \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$$

3. 分母 $\sum_{j=1}^{K} \pi_j \mathcal{N}(x_i \mid \mu_j, \Sigma_j)$：把所有声部的吸引力加起来，相当于这团音符的“总混音响度”。

$$\sum_{j=1}^{K} \pi_j \mathcal{N}(x_i \mid \mu_j, \Sigma_j)$$

4. 归一化：分子除以分母后，所有责任度之和等于 1，变成一套合法的概率。

$$\sum_{k=1}^{K} \gamma_{ik} = 1$$

调音师在分轨时不会武断“一刀切”，而是给每个音符写一张“声部来源概率表”。GMM 的数据点也一样：它属于每个簇的概率写满一整行，加起来恰好 100%。这张概率表，正是 EM 算法下一步腾挪的资本。

EM算法求解

调音师怎么从零开始分离声部？他两手交替：先用当前的乐器参数猜每个音符属于哪个声部，再用猜好的归属重新估计每件乐器的音高、音色和音量，然后拿着新参数再猜一遍……循环往复，混音越来越清晰。这个交替过程就是 EM 算法——期望最大化。

E 步（Expectation）：固定参数，用软聚类公式算出每个点的责任度 γ。M 步（Maximization）：固定责任度，用加权平均更新均值 μ、协方差 Σ 和权重 π。每一步都保证“对数似然”只升不降。为什么不能直接求最优解？因为对数似然里 log 套着求和，直接求导没有闭式解，只能像调音师一样一步步逼近。

$$\mu_k^{\text{new}} = \frac{\sum_{i=1}^{N}\gamma_{ik} x_i}{\sum_{i=1}^{N}\gamma_{ik}}, \qquad \pi_k^{\text{new}} = \frac{1}{N}\sum_{i=1}^{N}\gamma_{ik}, \qquad \Sigma_k^{\text{new}} = \frac{\sum_{i=1}^{N}\gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_{i=1}^{N}\gamma_{ik}}$$

1. 新均值 $\mu_k^{\text{new}}$：以责任度为权重的加权平均——每个点“投票”给这个声部的比例，决定声部中心往哪挪。

$$\mu_k^{\text{new}} = \frac{\sum_{i=1}^{N}\gamma_{ik} x_i}{\sum_{i=1}^{N}\gamma_{ik}}$$

2. 分子 $\sum_{i=1}^{N}\gamma_{ik} x_i$：所有数据点按责任度加权的总和，责任大的点影响更大。

$$\sum_{i=1}^{N}\gamma_{ik} x_i$$

3. 分母 $\sum_{i=1}^{N}\gamma_{ik}$：这个声部收到的总责任度，相当于“有多少音符主要归它管”。

$$\sum_{i=1}^{N}\gamma_{ik}$$

4. 新权重 $\pi_k^{\text{new}}$：声部总责任占全体数据的比例，等于重新校准音量推子。

$$\pi_k^{\text{new}} = \frac{1}{N}\sum_{i=1}^{N}\gamma_{ik}$$

5. 新协方差 $\Sigma_k^{\text{new}}$：按责任度加权统计散开的程度，决定椭圆的长短轴与朝向。

$$\Sigma_k^{\text{new}} = \frac{\sum_{i=1}^{N}\gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_{i=1}^{N}\gamma_{ik}}$$

注意 M 步里的责任度全部来自 E 步，两步骤相互依赖，这就是“交替迭代”的精髓。EM 保证对数似然单调不降，但不保证全局最优——调音师如果一开始把长笛听成低音提琴，后面很可能一直错下去。因此实践中往往多次随机初始化，取似然最高的结果。

典型应用场景

调音师学成出师后，接了三种活儿：把整段录音拆成各声部（软聚类）；给客户描述这场音乐会的“音色分布特征”（密度估计）；以及——听到一声不属于任何乐器的刺耳杂音，立刻警觉“器材出问题了”（异常检测）。一套 GMM 训练完毕，这三个能力自动到位。

软聚类：给每个新样本计算责任度 γ，得到概率归属。密度估计：直接计算 p(x)，知道任意位置的数据有多常见。异常检测：正常样本密度高，异常样本落在低密度角落，设定阈值即可判断。工业界大量使用这套组合拳，尤其是在传感器数据上。

$$\text{score}(x_i) = \sum_{k=1}^{K}\pi_k \, \mathcal{N}(x_i \mid \mu_k, \Sigma_k), \qquad x_i \text{ 异常} \iff \text{score}(x_i) < \varepsilon$$

1. 密度得分 $\text{score}(x_i)$：该点在混音模型里的“响度”，正常点响度高，异常点几乎无声。

$$\text{score}(x_i) = \sum_{k=1}^{K}\pi_k \, \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$$

2. 单成分贡献 $\pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$：每个声部都对这团声音贡献一份密度，加权求和即综合得分。

$$\pi_k \, \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$$

3. 阈值判断 $\varepsilon$：提前设定一个最小合理密度，低于它就触发告警。

$$\text{score}(x_i) < \varepsilon$$

一个正常音符在声学空间里密度很高，一声爆音则孤零零落在低密度区。GMM 的密度公式就是调音师手里的“分贝计”——读数越低，越可疑。这种思路在推荐、风控、运维监控里被反复验证，性价比极高。

与K-means对比

把聚类比作整理一堆混在一起的乐器声：K-means 像一位图省事的实习编辑，拿到每个音符只看哪个声部中心离它最近，一刀切归给那一个，绝不犹豫；GMM 的调音师则更细腻，他会记录每个音符对各声部的归属概率——同一个音符，可以同时写进小提琴和长笛两个分轨。

从数学上看，K-means 等价于“各成分方差相同、各向同性、权重相等”的特例；GMM 则允许每个簇有自己独立的协方差矩阵，形状可以是任意方向的椭圆。把硬分配的 argmin 换成软分配的后验概率 γ，K-means 就被“软化”成了 GMM。

$$c_i = \arg\min_{k}\|x_i - \mu_k\|^2 \qquad \text{vs} \qquad \gamma_{ik} = \frac{\pi_k\mathcal{N}(x_i\mid\mu_k,\Sigma_k)}{\sum_{j=1}^{K}\pi_j\mathcal{N}(x_i\mid\mu_j,\Sigma_j)}$$

1. 硬分配 $c_i$：K-means 里每个点只认一个最近的质心，输出单个标签。

$$c_i = \arg\min_{k}\|x_i - \mu_k\|^2$$

2. 欧氏距离 $\|x_i - \mu_k\|^2$：各方向等权计算距离，因此 K-means 只能画“圆”。

$$\lVert x_i - \mu_k \rVert^2$$

3. 软分配 $\gamma_{ik}$：GMM 里每个点是一条概率向量，对每个簇都有隶属度。

$$\gamma_{ik} = \frac{\pi_k\mathcal{N}(x_i\mid\mu_k,\Sigma_k)}{\sum_{j=1}^{K}\pi_j\mathcal{N}(x_i\mid\mu_j,\Sigma_j)}$$

实习编辑的快刀适合时间紧的粗剪，调音师的概率账本适合讲究的混音工程。选型经验很简单：数据量大、簇近似球形、只要硬标签，用 K-means；需要概率归属、簇呈椭圆或相互重叠、还要密度估计，就上 GMM。工具没有好坏，关键看你要的是“一个答案”还是“一组带着可信度的答案”。

AI 应用场景

Python 演示

import numpy as np  def gauss_pdf(x, mu, sigma):    """一元高斯概率密度: 钟形山的高度公式."""    coef = 1.0 / np.sqrt(2.0 * np.pi * sigma * sigma)    return coef * np.exp(-(x - mu) ** 2 / (2.0 * sigma * sigma))  def em_gmm_1d(x, K=3, max_iter=120):    x = np.asarray(x, dtype=float)    n = x.size    # 初始化: 沿数据分位点铺开 K 个声部中心    mu = np.quantile(x, np.linspace(0.05, 0.95, K))    sigma = np.full(K, x.std())    pi = np.full(K, 1.0 / K)    for _ in range(max_iter):        # E步: 计算责任度矩阵 resp (n x K)        resp = np.array([pi[k] * gauss_pdf(x, mu[k], sigma[k]) for k in range(K)]).T        resp /= resp.sum(axis=1, keepdims=True)        # M步: 按责任度加权更新参数        nk = resp.sum(axis=0)        mu = (resp.T @ x) / nk        dist2 = np.array([(x - mu[k]) ** 2 for k in range(K)]).T        sigma = np.sqrt((resp * dist2).sum(axis=0) / nk)        pi = nk / n    return mu, sigma, pi, resp  if __name__ == "__main__":    rng = np.random.default_rng(42)    # 三个真实声部: (均值, 标准差, 权重)    parts = [(0.0, 0.6, 0.4), (4.0, 0.8, 0.3), (8.0, 0.5, 0.3)]    x = np.concatenate([rng.normal(m, s, int(500 * w)) for m, s, w in parts])    mu, sigma, pi, resp = em_gmm_1d(x, K=3)    print("估计的均值:", np.round(mu, 2))    print("估计的标准差:", np.round(sigma, 2))    print("估计的权重:", np.round(pi, 2))    print("软聚类示例(前3个点的责任度):", np.round(resp[:3], 2))

1. coef = 1.0 / np.sqrt(2.0 * np.pi * sigma * sigma)钟形山的归一化系数，保证密度积分面积为 1。
2. resp = np.array([pi[k] * gauss_pdf(x, mu[k], sigma[k]) for k in range(K)]).TE 步：算出每个点对每个声部的原始“吸引力”。
3. resp /= resp.sum(axis=1, keepdims=True)按行归一化，让每个点的责任度加起来等于 1。
4. mu = (resp.T @ x) / nkM 步：责任度为权重的加权平均，更新声部中心。
5. sigma = np.sqrt((resp * dist2).sum(axis=0) / nk)M 步：按责任度加权统计散开程度，更新音色宽度。
6. pi = nk / nM 步：声部总责任占全体比例，更新音量推子。

常见误区

面试题