当前位置：首页 > news >正文

MATLAB实现GMM高斯混合模型聚类：从原理到代码实战

news 2026/9/1 6:35:56

简介：基于MATLAB编写的高斯混合模型（GMM）实现代码包，面向图像处理与视频分析方向的开发者，重点解决连续图像序列中的背景减除与前景检测问题。代码涵盖GMM建模流程，并通过EM算法完成参数估计，结构上分为模型初始化、密度计算、参数迭代更新与结果绘图等模块，便于读者对照原理逐段理解，也可直接用于监控场景下的运动目标提取实验。资源共7个文件，以5个M脚本为主，另含1个EPS示意图和1个MAT数据文件，压缩包整体约40KB，轻量精炼，适合快速阅读与二次开发。已有371人学习下载，可作为学习GMM原理、MATLAB算法实现及背景建模入门的一条高效路径。

1. GMM到底解决了什么问题

我先说一个自己的经历。早几年做用户行为聚类，一堆K-Means跑下去，聚类结果怎么看怎么别扭：用户群之间边界模模糊糊，有些点明明离两个簇中心都差不远，硬被切到一边，结果就是每个簇里都混着不少"异类"。后来换成了GMM（Gaussian Mixture Model，高斯混合模型），效果一下子顺了眼。这就是我想好好聊聊这套MATLAB写的GMM代码的起因。

GMM和K-Means最大的区别在于：K-Means是硬分配，一个点必须属于某一个簇；GMM是软分配，一个点可以以不同概率属于多个簇。拿现实打比方，K-Means像班主任硬性分组，每个学生只能进一个兴趣小组；GMM则更像同时给多个社团投票，按比例分配你的时间和精力。很多数据本身就不是一刀切能切开的，比如收入分层、图像像素分类、语音特征聚类，天然就带着"重叠区"，这时候GMM就明显比K-Means靠谱。

这套MATLAB代码适合谁来参考？正在做聚类分析、密度估计、数据建模的研究生和工程师，或者想把手写EM算法跑通的机器学习初学者。它能帮你不依赖现成工具箱，从底层理解GMM的运行机制，同时也提供了直接可用的高效版本，方便快速落地到具体任务里。

2. 数学原理不绕弯子

2.1 高斯混合模型是个什么结构

所谓混合模型，就是把若干个高斯分布按一定权重叠加在一起。二维高斯分布就像一座山，有山顶（均值）、有胖瘦（协方差矩阵）；多个高斯分布加权求和，就得到一片连绵起伏的"山脉"。每个样本来自哪座山是未知的，这就是隐变量。

数学上写成：

p(x) = Σₖ πₖ · N(x | μₖ, Σₖ)

其中πₖ是第k个分量的权重，满足Σπₖ=1，N(x | μₖ, Σₖ)是均值为μₖ、协方差为Σₖ的高斯分布。GMM要做的事，就是给定一堆样本X，反推出最可能的πₖ、μₖ、Σₖ。

2.2 EM算法简化拆解

直接对似然函数求导没法解出解析解，因为每个样本属于哪个分量是未知的。EM算法绕开了这个难题，分两步交替迭代：

E步（Expectation）：在当前参数下，计算每个样本属于每个分量的后验概率，也就是"责任度"γₙₖ。这一步相当于先猜：按现在的参数，这些点大概来自哪些山？

M步（Maximization）：用责任度做权重，重新估计参数。权重πₖ是责任度的平均，μₖ是责任度加权的样本均值，Σₖ是责任度加权的协方差。这一步相当于根据猜测结果，重新把每座山的位置和形状修一遍。

E步和M步交替迭代，每轮都会让对数似然增加，直到收敛。收敛后你就得到了一组参数，可以直接用来做聚类（取最大责任度对应的分量）或者做密度估计（代入p(x)计算任意点的概率密度）。

3. 完整MATLAB代码实现

3.1 主函数结构与初始化策略

先给出完整的主函数。这段代码我写成通用形式，输入数据矩阵X（N×D），聚类数K，输出各个分量的参数和聚类标签：

function [label, mu, Sigma, pi, logL] = my_gmm(X, K, maxIter) % 手写GMM(高斯混合模型)聚类 % 输入: % X: N×D 数据矩阵 % K: 聚类数 % maxIter: 最大迭代次数, 默认200 % 输出: % label: N×1 聚类标签 % mu: K×D 均值 % Sigma: D×D×K 协方差 % pi: K×1 权重 % logL: 每轮对数似然值 if nargin < 3 maxIter = 200; end [N, D] = size(X); % 用K-Means先粗聚类, 得到初始参数 [~, initLabel] = kmeans(X, K, 'MaxIter', 50, 'Replicates', 3); % 初始化均值 mu = zeros(K, D); Sigma = zeros(D, D, K); pi = zeros(K, 1); for k = 1:K idx = (initLabel == k); if sum(idx) > 1 mu(k, :) = mean(X(idx, :), 1); Sigma(:, :, k) = cov(X(idx, :)) + 1e-6 * eye(D); else mu(k, :) = X(randi(N), :); Sigma(:, :, k) = eye(D); end pi(k) = sum(idx) / N; end logL = zeros(maxIter, 1); % 存储责任度 gamma = zeros(N, K); for iter = 1:maxIter % ---- E步: 计算责任度 ---- for k = 1:K gamma(:, k) = pi(k) * mvnpdf(X, mu(k, :), Sigma(:, :, k)); end % 归一化 gamma_sum = sum(gamma, 2); gamma = gamma ./ repmat(gamma_sum, 1, K); % ---- 计算对数似然 ---- logL(iter) = sum(log(gamma_sum + eps)); % ---- M步: 更新参数 ---- Nk = sum(gamma, 1); for k = 1:K if Nk(k) < 1e-6 continue; % 避免分量消失 end mu(k, :) = (gamma(:, k)' * X) / Nk(k); Xc = X - repmat(mu(k, :), N, 1); Sigma(:, :, k) = (Xc .* repmat(gamma(:, k), 1, D))' * Xc / Nk(k); Sigma(:, :, k) = Sigma(:, :, k) + 1e-6 * eye(D); pi(k) = Nk(k) / N; end % 检查收敛 if iter > 1 && abs(logL(iter) - logL(iter-1)) < 1e-6 logL = logL(1:iter); break; end end % 分配标签: 取责任度最大的分量 [~, label] = max(gamma, [], 2); end

kmeans初始化这一步非常关键。随机初始化对GMM不友好，容易掉进局部最优，而K-Means先跑一遍，能得到一个相对合理的初始位置，后续EM收敛又快又稳。我给协方差矩阵加了1e-6的对角扰动，这是为了防止出现奇异矩阵导致mvnpdf报错。

3.2 注意防止数值下溢

E步里有个隐形坑：当数据维度较高或者分量距离较远时，pi(k) * mvnpdf(...)算出来的值可能极小，甚至小到MATLAB都分辨不出来，直接变0。一旦某一行的所有分量都是0，归一化就会得到NaN，程序直接罢工。eps的添加能兜住求对数时的无穷大问题，但更稳妥的办法是每轮都检查gamma_sum里有没有0，有就做一次平移处理。实际项目中我建议用log-sum-exp技巧改写E步，能从根本上避免这个问题。

3.3 生成测试数据观察聚类效果

跑通逻辑最直接的办法是生成一份已知分布的数据来验证。比如生成三个有部分重叠的高斯簇，然后看聚类结果准不准：

% 生成立即可视化测试数据 rng(42); N1 = 300; N2 = 400; N3 = 300; X1 = mvnrnd([0 0], [1 0.3; 0.3 1.2], N1); X2 = mvnrnd([4 3], [1.5 0; 0 1], N2); X3 = mvnrnd([2 6], [0.8 0.2; 0.2 1.1], N3); X = [X1; X2; X3]; [label, mu, Sigma, pi, logL] = my_gmm(X, 3, 300); figure; gscatter(X(:,1), X(:,2), label, 'rgb', 'o', 8); hold on; plot(mu(:,1), mu(:,2), 'kx', 'MarkerSize', 12, 'LineWidth', 2); title('手写GMM聚类结果'); legend('簇1', '簇2', '簇3', '聚类中心'); grid on;

跑完之后你会看到，三个簇之间原本交叠的点被合理地分到了概率更大的一侧，不像K-Means那样直接画一条硬边界。再看logL曲线，基本在10到20轮内就稳定下来，说明收敛速度相当理想。

4. 用MATLAB官方工具箱函数做对比验证

4.1 fitgmdist一行命令搞定

手写代码是为了搞懂原理，真正工程落地时直接调用fitgmdist更高效。MATLAB统计工具箱自带的高斯混合模型拟合函数，用法极其简洁：

% 官方工具箱版本 gm = fitgmdist(X, 3, 'RegularizationValue', 1e-5); label = cluster(gm, X); % 查看参数 gm.mu gm.Sigma gm.ComponentProportion

fitgmdist底层实现比手写版本稳健得多，它自带多种初始化策略、正则化处理和边界约束，还会自动处理分量退化问题。对大多数应用来说，直接调这个函数就够了。

4.2 手写版与工具箱版的取舍

对比项	手写版	fitgmdist
代码量	约60行	1行
可定制性	高，可改任何步骤	中，只能调选项参数
数值稳定性	需自己加保护	内置正则化
学习价值	极高	一般
运行速度	中	快（C底层）
依赖工具箱	仅基础MATLAB	Statistics Toolbox

如果只是做数据分析、聚类实验，直接用fitgmdist；如果是为了课程设计、面试准备或二次开发定制算法，建议认真读一遍手写实现。我自己做项目时，会先用工具箱版本跑通流程，再用定制版处理特殊需求，比如给某些分量固定均值的先验约束。

5. K值怎么选，协方差结构怎么定

5.1 用BIC/ AIC 选择聚类数

GMM最头疼的就是K值的确定。K-Means用肘部法则，GMM则更常看信息准则。思路是在不同K下分别拟合模型，看哪个K能在"拟合度"和"复杂度"之间取得平衡。BIC对复杂度惩罚更狠，通常选择的K更小更简洁：

K_candidates = 1:8; BIC_values = zeros(size(K_candidates)); AIC_values = zeros(size(K_candidates)); for i = 1:length(K_candidates) gm = fitgmdist(X, K_candidates(i), 'RegularizationValue', 1e-5); AIC_values(i) = gm.AIC; BIC_values(i) = gm.BIC; end figure; plot(K_candidates, BIC_values, 'o-', 'LineWidth', 1.5); hold on; plot(K_candidates, AIC_values, 's--', 'LineWidth', 1.5); xlabel('K值'); ylabel('信息准则值'); legend('BIC', 'AIC'); grid on; % 找BIC最小值对应的K [~, bestK] = min(BIC_values); fprintf('BIC最优K = %d\n', bestK);

选K的时候有个细节：BIC曲线往往不会只有一个孤零零的最小值，而是先骤降再缓慢爬升，此时选"拐点"比选"全局最小"更实用。因为全局最小K有时会过拟合噪声，拐点K则保留了更强的泛化能力。

5.2 协方差结构的三种选择

fitgmdist里的CovarianceType参数值得专门提一下：

full（全协方差）：每个分量有自己完整的协方差矩阵，能拟合不同形状的簇，但参数多、容易过拟合，需要足够样本量。
diag（对角协方差）：每个分量各维度独立，参数大幅减少，拟合速度更快，适合高维数据或维度间相关性弱的场景。
tied（共享协方差）：所有分量共用同一个协方差矩阵，相当于大家形状一样、大小差不多，只是位置不同。适合簇形状高度一致的场景，参数最少。

实际用的时候，先跑一遍diag看效果，如果聚类结果明显不合理，再升级为full。直接上full在高维数据上极易出状况，训练慢还不稳定。

6. 高频踩坑与排查技巧

6.1 协方差矩阵奇异

这是GMM最经典的问题。当某个分量的样本太少或维度太高时，协方差矩阵可能变成奇异矩阵，mvnpdf直接报错。

现象：报错"SIGMA must be symmetric positive definite"。

原因：一个分量的责任度几乎全为0，导致M步里Nk(k)极小，协方差估计失真。

排查方案：打印每轮的Nk，看哪些分量快要"饿死"；检查数据是否有完全共线的维度（比如某列是另一列的2倍）；把样本量除以维度数，如果比值小于5，建议降维或改用diag协方差。

6.2 代码能跑但聚类效果差

有时候程序不报错，但聚类结果明显不对，这种情况多数是初始化没做好。

经验做法：固定随机种子；多次重置初始值跑多个起点，选对数似然最高的结果；手写版里把K-Means的Replicates从3提高到10，能显著提升稳定性。我自己踩过最坑的一次是数据里有两个簇中心几乎重合，K-Means初始化时把两个中心都扔到了同一团点附近，结果EM迭代到最后两个分量完全重叠，白白浪费了一整轮实验。加大Replicates之后基本不会再遇到这种问题。