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

简介基于MATLAB编写的高斯混合模型GMM实现代码包面向图像处理与视频分析方向的开发者重点解决连续图像序列中的背景减除与前景检测问题。代码涵盖GMM建模流程并通过EM算法完成参数估计结构上分为模型初始化、密度计算、参数迭代更新与结果绘图等模块便于读者对照原理逐段理解也可直接用于监控场景下的运动目标提取实验。资源共7个文件以5个M脚本为主另含1个EPS示意图和1个MAT数据文件压缩包整体约40KB轻量精炼适合快速阅读与二次开发。已有371人学习下载可作为学习GMM原理、MATLAB算法实现及背景建模入门的一条高效路径。1. GMM到底解决了什么问题我先说一个自己的经历。早几年做用户行为聚类一堆K-Means跑下去聚类结果怎么看怎么别扭用户群之间边界模模糊糊有些点明明离两个簇中心都差不远硬被切到一边结果就是每个簇里都混着不少异类。后来换成了GMMGaussian 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个分量的权重满足Σπₖ1N(x | μₖ, Σₖ)是均值为μₖ、协方差为Σₖ的高斯分布。GMM要做的事就是给定一堆样本X反推出最可能的πₖ、μₖ、Σₖ。2.2 EM算法简化拆解直接对似然函数求导没法解出解析解因为每个样本属于哪个分量是未知的。EM算法绕开了这个难题分两步交替迭代E步Expectation在当前参数下计算每个样本属于每个分量的后验概率也就是责任度γₙₖ。这一步相当于先猜按现在的参数这些点大概来自哪些山M步Maximization用责任度做权重重新估计参数。权重πₖ是责任度的平均μₖ是责任度加权的样本均值Σₖ是责任度加权的协方差。这一步相当于根据猜测结果重新把每座山的位置和形状修一遍。E步和M步交替迭代每轮都会让对数似然增加直到收敛。收敛后你就得到了一组参数可以直接用来做聚类取最大责任度对应的分量或者做密度估计代入p(x)计算任意点的概率密度。3. 完整MATLAB代码实现3.1 主函数结构与初始化策略先给出完整的主函数。这段代码我写成通用形式输入数据矩阵XN×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); endkmeans初始化这一步非常关键。随机初始化对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.ComponentProportionfitgmdist底层实现比手写版本稳健得多它自带多种初始化策略、正则化处理和边界约束还会自动处理分量退化问题。对大多数应用来说直接调这个函数就够了。4.2 手写版与工具箱版的取舍对比项手写版fitgmdist代码量约60行1行可定制性高可改任何步骤中只能调选项参数数值稳定性需自己加保护内置正则化学习价值极高一般运行速度中快C底层依赖工具箱仅基础MATLABStatistics 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之后基本不会再遇到这种问题。6.3 对数似然震荡不收敛现象logL曲线不是单调上升而是上下跳动。原因数据量太小、异常值太多或者E步出现了NaN没有报错。排查方案在E步后加一行检查any(isnan(gamma(:)))有就停下来打印当前参数对数据进行标准化让每个维度均值0方差1数值稳定性会好很多检查代码里是否用了inv而不是\协方差求逆建议用pinv或\。6.4 分量消失问题EM迭代中某个分量的权重趋于0对应的高斯分布被挤掉了。这在数据本身只有K-1个簇时会自然发生不一定算bug但如果你确认数据应该有K个簇就要检查初始化是不是把某个中心放得太远。7. 三个水调输出前避坑我在项目里实际跑过不少GMM相关的任务最后再分享几个让代码更抗造的小改动。第一个改动是数据标准化。输入到GMM之前先对每一列做z-score标准化减均值除标准差。GMM对量纲极其敏感某个维度的数值范围稍微大一点就会主导距离计算让其他维度的贡献几乎为零。标准化的代价是聚类结果不好直接解释但如果你只关心分组结构标准化绝对值得。第二个改动是尽量用cluster(gm, X)而不是手动max(gamma, [], 2)拿标签。尽管两者理论上等价工具箱的cluster额外包含了处理NaN和边缘情况的安全逻辑线上环境更可靠。第三个改动是画图时把后验概率可视化出来。不要只画硬标签散点图试着用透明度或颜色深浅表示每个点属于目标簇的最大概率。你经常会发现那些边界区的点概率只有0.4、0.5这时候硬聚类结果带来的好像分得很清的错觉就会消失也会提醒你数据本身可能确实没有清晰的分群结构。最后有一个忠告GMM不是万能聚类器。如果数据呈现明显的非凸形状比如环形、螺旋形GMM的椭圆簇假设会让它无能为力这时候考虑DBSCAN或其他密度聚类方法更合适。但在凸簇、部分重叠、密度不均的常规场景里GMM仍然是综合表现最稳的模型之一。手写一遍你能真正理解它擅长什么也才知道它怕什么。本文还有配套的精品资源点击获取