MATLAB手写GMM：EM算法实现、避坑指南与参数调优

简介：面向图像处理与计算机视觉学习者的MATLAB高斯混合模型（GMM）实现代码，重点演示如何用GMM完成视频背景减除与前景目标提取。代码采用期望最大化（EM）算法，完整覆盖高斯分量初始化（k均值聚类）、概率密度计算、参数迭代更新与前景判定等关键步骤，便于读者对照教材理解原理，也可直接嵌入自身检测项目验证效果。压缩包共7个文件，其中5个.m脚本构成核心算法与绘图模块，1个.mat文件存放实验输入数据，1个.eps文件输出结果图，整体仅40KB，结构清晰、注释精简，适合快速阅读和二次开发。该资源已有373人学习下载，对想快速上手GMM背景建模的读者来说，是一份体量小巧且可直接运行的实战样例，能有效节省从头搭建算法框架的时间，同时帮助深入理解GMM在图像序列处理中的工程实现细节。

1. matlab写的GMM代码：别急着调工具箱，从EM算法开始反而更快

很多人在MATLAB里做GMM第一反应就是调fitgmdist，但等你想改一点东西——比如让某个分量的协方差固定、加一个权重约束、或者在新数据进来时做增量更新——工具箱的黑匣子就卡住你了。我做过的模拟项目X里，早期就直接用fitgmdist出聚类结果，直到有一天一个分量坍缩到单个样本点上，我才把matlab写的GMM代码整个扒开重写了一遍。实际上，GMM的本质就是EM算法在混合高斯上的具体展开，自写一遍的代码量并不大，核心E步和M步加起来不过几十行。对想搞懂模型行为、想定制算法的从业者来说，自己写反而比调包更快，而且后面的所有参数调节都不是瞎猜的。

这篇就把我常用的那套matlab写的GMM代码拆开讲清楚：从数学骨架、向量化实现、初始化策略，到收敛判断和避坑。新手能照着跑通一个小例程，熟手能在参数和边界上找到参考。最后再给一个BIC选分量数的收尾技巧。

2. GMM的数学骨架与MATLAB中的向量化建模

2.1 高斯混合模型的核心假设与分量设计

GMM假设你的数据是从若干个高斯分量里混合生成的。对一个D维样本x，概率密度写成：

p(x) = sum(k=1..K) pi_k * N(x | mu_k, Sigma_k)

其中pi_k是混合系数，满足pi_k >= 0且sum(pi_k)=1。每个分量是一个D维高斯分布，由均值mu_k和协方差Sigma_k决定。你要估计的东西就是这三组参数：mu、Sigma、pi。

参数个数比较容易算：每个分量的均值有D个自由度，协方差矩阵有D*(D+1)/2个自由度，混合系数因为和为1约束，只有K-1个自由度。总参数个数是KD + KD*(D+1)/2 + (K-1)，这个数字后面算BIC要用，现在先记住。

K本身不是EM跑出来的，是你要在运行前定的超参数。常见做法是跑K=1到K=8，算BIC或AIC来选。注意，如果你只用极大似然不加惩罚，K越大似然越高，最后会过拟合到每个点一个分量，所以模型选择这一步绕不开。

2.2 从似然函数到EM迭代：用矩阵运算改写两个关键步骤

GMM的极大似然没有解析解，原因是log里面套着求和，直接求导等于零解不出来。EM的思路是引入隐变量——每个样本属于哪个分量——然后迭代优化一个下界。E步算后验概率，M步用这些后验概率做加权最大似然更新参数。

对N个样本、K个分量，E步的核心是算responsibility矩阵（通常叫R或gamma），形状是N×K，每个元素R(n,k)表示第n个样本属于第k个分量的概率：

R(n,k) = pi_k * N(x_n | mu_k, Sigma_k) / sum_j pi_j * N(x_n | mu_j, Sigma_j)

M步更新公式：

mu_k = sum_n R(n,k) * x_n / sum_n R(n,k)

Sigma_k = sum_n R(n,k) * (x_n - mu_k)(x_n - mu_k)^T / sum_n R(n,k)

pi_k = sum_n R(n,k) / N

在MATLAB里用for循环写这三个更新其实也能跑，但N到几万、K到十来个的时候，循环会慢得让人想放弃。我会把所有操作改成矩阵运算：样本矩阵X是N×D，均值矩阵mu是K×D，协方差用三维数组Sigma(:,:,k)存。E步里最费时间的部分是马氏距离计算，不要每个分量循环算完再拼接，直接向量化后拼成N×K矩阵。

2.3 MATLAB向量化实现E步与M步的核心代码

先给一个不依赖统计工具箱的完整函数骨架，这版是我现在改项目时最常用的底子：

function [mu, Sigma, p, loglik] = gmm_em_custom(X, K, maxIter, tol, reg) % 自写GMM-EM训练函数，不依赖统计工具箱 % 输入: % X - NxD 数据矩阵，N为样本数，D为维度 % K - 高斯分量个数 % maxIter - 最大迭代次数，默认200 % tol - 对数似然变化量的收敛阈值，默认1e-6 % reg - 协方差对角正则项，防止奇异，默认1e-6 % 输出: % mu - KxD 均值矩阵 % Sigma - DxDxK 协方差数组 % p - Kx1 混合系数 % loglik - 最终对数似然值 [N, D] = size(X); % ---------- 初始化 ---------- rng(42); % 固定随机种子便于复现 idx = randperm(N, K); mu = X(idx, :); Sigma = repmat(cov(X), 1, 1, K); % 初始协方差用全局协方差 p = ones(K, 1) / K; prev_loglik = -inf; for iter = 1:maxIter % ---------- E步 ---------- resp_raw = zeros(N, K); for k = 1:K diff = X - mu(k, :); % 用 \ 解线性方程，避免显式求逆 mahal = sum((diff / Sigma(:, :, k)) .* diff, 2); detSigma = det(Sigma(:, :, k)); resp_raw(:, k) = p(k) * exp(-0.5 * mahal) / sqrt((2*pi)^D * detSigma); end resp = resp_raw ./ sum(resp_raw, 2); % 归一化得到后验概率 % ---------- M步 ---------- Nk = sum(resp, 1); % 每个分量的有效点数，1xK mu = (resp' * X) ./ Nk'; for k = 1:K diff = X - mu(k, :); Sigma(:, :, k) = (diff' * (resp(:, k) .* diff)) / Nk(k); % 加正则项，防止奇异协方差 Sigma(:, :, k) = Sigma(:, :, k) + reg * eye(D); end p = Nk' / N; % ---------- 对数似然与收敛判断 ---------- loglik_new = sum(log(sum(resp_raw, 2))); % 注意用未归一化的resp_raw if abs(loglik_new - prev_loglik) < tol loglik = loglik_new; break; end prev_loglik = loglik_new; end loglik = prev_loglik; end

这段代码里有几个地方容易踩坑，我说一下我的习惯。首先，E步里resp_raw是每个样本的原始加权高斯密度值，resp是归一化后的后验概率。计算对数似然时，一定用resp_raw的行和取log再求和，因为resp的行和恒等于1，拿它算loglik就永远是0，这是个很隐蔽的错。其次，马氏距离我写成diff / Sigma再乘diff，也就是diff * inv(Sigma) * diff'的向量化形式，但用的是矩阵右除而不是inv，数值稳定性好一些。

协方差初始化为全局cov(X)，而不是单位矩阵，能让第一轮迭代的loglik不至于离谱，收敛也快很多。复制到每个分量时用了repmat扩展成D×D×K的三维数组，后续更新Sigma时直接对第k个平面赋值。reg参数我默认给1e-6，但遇到数据维度高或样本分布窄时，要调到1e-4甚至1e-2，下面避坑章里细说。

3. 搭一个可跑的MATLAB GMM训练脚本：从模拟数据到聚类结果

3.1 样本初始化：选点、协方差、混合系数的三个讲究

初始化决定了EM收敛到哪个局部最优，这是我调试时最常调的一块。三个选择都有讲究：

第一个是均值mu的初始化。最土的办法是随机选K个样本点当初始均值，我一般会用randperm选不重复的索引。但纯随机在分量靠近时容易初始化到同一个簇里，所以更稳的做法是先跑一次kmeans，把kmeans的聚类中心拿来当EM的初始均值。

第二个是协方差Sigma的初始化。不要用单位矩阵，也不要每个分量随机生成，不然前几轮迭代的密度值会非常扭曲。我习惯用全样本的协方差矩阵去赋值给每个分量，只加一点小扰动。这样可以保证初始时每个高斯分量都覆盖整个数据区域，EM在早期会把分量慢慢拉开，而不是一上来就互相挤压。

第三个是混合系数p的初始化。没有特殊先验时，均匀分布就是最好的起点，即p = ones(K,1)/K。注意M步更新后p是Nk/N，加起来一定为1，不需要额外归一化。

% 初始化示例：kmeans版本 rng(42); % 先跑一次kmeans v，用聚类中心初始化mu [idx, C] = kmeans(X, K, 'MaxIter', 100); mu = C; Sigma = repmat(cov(X), 1, 1, K); p = ones(K, 1) / K;

这里用了kmeans需要统计工具箱，如果你也想完全不依赖工具箱，那就用纯随机选点那版。两种效果差距不算大，但kmeans初始化能显著减少EM的迭代次数。我在模拟项目X里做过对比：纯随机初始化平均要跑80轮收敛，kmeans初始化平均30轮就到了。

3.2 主循环：迭代更新mu、Sigma与混合系数的MATLAB写法

主循环就是第2章的代码骨架在跑，但要加一个保险：记录每一轮的对数似然变化。这个序列非常有用，既能判断是否收敛，也能排查发散问题。我一般在函数外再包一层，把loglik序列存下来画图：

loglik_history = zeros(maxIter, 1); for iter = 1:maxIter % E步、M步、收敛判断代码同前 loglik_history(iter) = loglik_new; if abs(loglik_new - prev_loglik) < tol loglik_history = loglik_history(1:iter); break; end end % 画对数似然曲线，排查发散或收敛过慢 figure; plot(loglik_history, 'o-'); xlabel('迭代次数'); ylabel('对数似然'); title('GMM-EM对数似然收敛曲线'); grid on;

loglik曲线正常形态是单调上升，然后在一段平台期趋于平稳。如果你看到曲线有下降段，说明代码里有bug，常见原因是用归一化后的resp算了似然，或者协方差更新后没有加正则项导致数值爆掉。如果你看到收敛特别快但聚类结果明显不对，那就回到初始化，换一组随机种子重新开始。

3.3 造一组合成数据测试：自己验证GMM代码是否可靠

没有真实数据时，我会用已知参数生成一组合成数据来验证代码。比如生成三个高斯分量，均值分别在(-3,-3)、(3,3)、(0,-2)，协方差各不相同，混入噪声点，然后跑自己的GMM函数：

% 生成合成数据 rng(1); N = 1000; % 三个分量，各占三分之一 X = [ mvnrnd([-3 -3], [0.5 0; 0 0.8], 300); mvnrnd([3 3], [1 0.3; 0.3 1], 400); mvnrnd([0 -2], [0.3 0; 0 0.3], 300) ]; % 调用自写GMM K = 3; [mu, Sigma, p, loglik] = gmm_em_custom(X, K, 200, 1e-6, 1e-5); % 可视化结果 figure; gscatter(X(:,1), X(:,2), cluster(X, mu, Sigma, p)); hold on; plot(mu(:,1), mu(:,2), 'kx', 'MarkerSize', 12, 'LineWidth', 2);

这里的cluster函数是拿E步算出来的后验概率做硬分类，也就是对每个样本取responsibility最大的分量编号：

function label = cluster(X, mu, Sigma, p) N = size(X, 1); K = size(mu, 1); resp = zeros(N, K); for k = 1:K diff = X - mu(k, :); mahal = sum((diff / Sigma(:, :, k)) .* diff, 2); detSigma = det(Sigma(:, :, k)); resp(:, k) = p(k) * exp(-0.5 * mahal) / sqrt((2*pi)^D * detSigma); end [~, label] = max(resp, [], 2); end

这套验证流程我每次都跑一遍，确认能恢复出原始分量的均值和大致协方差，才敢把代码用到新数据上。如果均值偏差超过1个标准差，基本可以断定是收敛到局部最优了，需要调整初始化或增加随机重试次数。

4. 与fitgmdist对比：什么时候自己写，什么时候用工具箱

4.1 自写代码与fitgmdist的功能差异

fitgmdist是MATLAB统计工具箱里封装好的函数，它解决的问题和自写代码完全一样，但实现里做了大量工程化处理。它默认用kmeans++初始化，还有正则化、协方差类型可选、并行计算、以及现成的BIC接口。

fitgmdist的关键参数如下表，这也是你评估要不要换掉自写代码的对照表：

用的时候一句话就够了：

gmm = fitgmdist(X, K, 'CovarianceType', 'full', ... 'RegularizationValue', 1e-5, 'Replicates', 10, ... 'Options', statset('Display', 'final', 'MaxIter', 200));

输出的是gmm对象，可以用gmm.mu、gmm.Sigma、gmm.ComponentProportion拿到参数，用gmm.posterior(X)算后验概率，用gmm.logL拿最终对数似然。它省掉的麻烦非常多，特别是Replicates重试机制，自己写要额外铺一坨代码。

4.2 我什么情况下坚持用自写代码

fitgmdist虽好，但它的对象结构和训练过程对某些场景不友好。我在三类场景下会坚持用自写GMM代码：

第一种是想让不同分量的协方差结构不同。fitgmdist要么所有分量都是full，要么都是diagonal，没法设置“分量1用对角、分量2用完整”。自写代码里Sigma三维数组想怎么赋值就怎么赋值，改一行就行。

第二种是增量学习。新数据一批一批来，不想把老数据全量存下来重新训练。自写代码可以拿旧的mu、Sigma、p当新EM的初始值，只迭代十几轮就收敛，而fitgmdist的Start参数虽然支持传自定义初始值，但内部实现还是全量数据训练。

第三种是嵌入到别的模型里。比如我在做某跨平台系统的在线聚类模块时，需要每5分钟拿新数据更新一次GMM，而且失败要能重启恢复。自写代码可以输出每一步的完整结构，甚至能把中间状态存成MAT文件，fitgmdist对象序列化不如直接存数组方便。

4.3 工具箱兜底，自写代码做验证的联合用法

我自己现在的习惯是两者都写：自写代码用于开发期调试，fitgmdist作为验证基准。每次改完EM的更新公式或初始化策略，都会跑一遍同一份数据做对照：

% fitgmdist作为基准 gmm_ref = fitgmdist(X, K, 'RegularizationValue', 1e-5, 'Replicates', 10); ref_loglik = gmm_ref.logL; % 自写代码 [mu, Sigma, p, my_loglik] = gmm_em_custom(X, K, 200, 1e-6, 1e-5); % 对比对数似然，误差在0.1%以内视为正常 fprintf('fitgmdist loglik: %.4f\n', ref_loglik); fprintf('custom loglik: %.4f\n', my_loglik);

如果自写代码的loglik比fitgmdist差太多，优先怀疑收敛判断或者初始化出问题，而不是算法公式出错。反过来讲，如果你只是要一个聚类结果，不在乎内部参数的行为，fitgmdist绝对够了。不要因为“自己写更酷”就重复造轮子，工程时间也是成本。

5. 避坑指南：自写GMM在MATLAB里翻车的5个高频问题

5.1 协方差矩阵奇异导致loglik变NaN

现象：迭代到某一步后，对数似然变成NaN，或者出现警告“Matrix is singular to working precision”。

原因：某个分量的有效样本点数Nk太小，或者数据本身有强线性相关性，导致更新出的协方差矩阵半正定退化。比如二维数据里两个特征高度相关，协方差矩阵的行列式趋近于零，exp里面的马氏距离除以接近0的数就爆掉了。

解决：M步更新完协方差后，立刻加正则项。我的做法是Sigma(:, :, k) = Sigma(:, :, k) + reg * eye(D)，reg取值从1e-6开始，如果还奇异就提高到1e-4。另外可以在每次更新后检查det(Sigma(:, :, k))是否小于某个阈值，比如1e-10，是的话强制重置该分量的协方差为全局协方差cov(X)。千万不要把这个检查省掉，高维数据下这个坑几乎是必踩的。

5.2 一个分量坍缩到单个样本点

现象：聚类结果里某个分量只包含一两个点，而且它的协方差矩阵变得特别小，对数似然反而还在上升。

原因：EM对离群点（outlier）很敏感。如果初始化时某个均值落在了远离数据主体的位置，或者某个离群点被单独分了一个分量，EM会不断把这个分量的方差压缩到极端小，因为它正好能解释那一个点的极大似然。

解决：一种做法是在M步更新后发现Nk < 5就重新初始化这个分量，但更常见的做法是限制协方差的下界，即正则项设置得偏高，让任何一个分量都无法收缩到单点尺度。另一种做法是换初始化，我常用kmeans初始化来避免这种情况。实际调参时，遇到分量坍缩的轮次，我第一件事不是改正则项，而是换随机种子重新跑，大部分情况下能绕过去。如果重试很多次都坍缩，那说明有一个离群点非常远，建议先把离群点剔除再做GMM。

5.3 不同随机种子结果差异大

现象：固定K不变，换一个随机种子，聚类结果整个变了，连分量个数看起来都不一样。

原因：EM是爬坡算法，收敛到哪个局部最优受初始化影响很大。GMM的目标函数非凸，不同的初始均值组合会收敛到不同的解。

解决：多次随机初始化然后保留loglik最大的那次，这是最笨但最可靠的方法。我一般跑10到20次，配一个for循环：

best_loglik = -inf; for trial = 1:20 [mu_tmp, Sigma_tmp, p_tmp, loglik_tmp] = gmm_em_custom(X, K, 200, 1e-6, 1e-5); if loglik_tmp > best_loglik best_loglik = loglik_tmp; mu = mu_tmp; Sigma = Sigma_tmp; p = p_tmp; end end

如果20次重试里最优解依然和fitgmdist差很远，那就检查你的数据里是否真的存在K个可分离的分量，或者K设大了，有可能其中一个分量需要覆盖一个极不规则的形状。

5.4 收敛判断过早停止

现象：loglik没跑几步就不变了，但聚类结果的轮廓明显不对，换大maxIter也没有改善。

原因：我早期只用绝对变化量abs(loglik_new - prev_loglik) < tol来判断收敛，当N很大时，每个样本的对数似然平均变化极小，但整体参数其实还有明显的调整空间。比如N=10万，tol=1e-4，其实可能只差几十的绝对变化，但这时参数还在漂。

解决：改用相对变化量判断：abs(loglik_new - prev_loglik) / abs(prev_loglik) < tol。或者同时跟踪参数变化，比如计算mu矩阵的Frobenius范数变化。我的经验是tol设为1e-6比较稳，maxIter设200到300，如果到了迭代上限还在缓慢上升，说明收敛标准太严，或者数据本身不太适合GMM。另外检查一下loglik曲线在收敛前是否有明显的平台期，如果有，说明算法卡住了，不是收敛。

5.5 二维可视化结果与预期不符

现象：画出的散点图和聚类边界比起来，肉眼可见地有一个簇被另一个簇“吃掉”，或者某个区域应该有一个分量但全被划到别处去了。

原因：GMM聚类的边界实际上是后验概率的等值线，不是简单的几何分割。如果你只画了硬分类的结果，两个分量重叠区域的样本被归到概率大的那边，看起来就会觉得“偏了”。另一个原因是数据本身的分量形状不是高斯，比如偏态分布或带长尾，高斯假设不成立时边界自然很怪。

解决：先画出后验概率等值线确认模型行为：

d = 0.1; [x1, x2] = meshgrid(min(X(:,1))-1:d:max(X(:,1))+1, min(X(:,2))-1:d:max(X(:,2))+1); Xgrid = [x1(:), x2(:)]; post = posterior_gmm(Xgrid, mu, Sigma, p); contour(x1, x2, reshape(post(:,1), size(x1)), [0.5 0.5], 'r');

posterior_gmm函数就是对每个网格点算归一化的responsibility。如果等值线显示两个分量其实分得很好，只是硬分类边界看着丑，那就没问题。如果等值线本身重叠严重，说明K设小了，或者数据不满足高斯假设，考虑换t混合模型或者加一个分量。这一条本质上是模型诊断，很多初学同学习惯跳过，但在我这里属于每次必做的一步。

6. 收尾：BIC定分量数与几个顺手就能用的扩展

自写GMM代码跑通之后，最实用的一件事就是用BIC确定分量数K。拟合完每个K，计算：

BIC = -2 * loglik + 参数个数 * log(N)

参数个数的公式前面提过：KD + KD*(D+1)/2 + (K-1)。MATLAB里按K循环跑一遍：

K_list = 1:8; bic_values = zeros(length(K_list), 1); loglik_values = zeros(length(K_list), 1); for i = 1:length(K_list) K = K_list(i); [mu, Sigma, p, loglik] = gmm_em_custom(X, K, 300, 1e-6, 1e-5); param_count = K*D + K*D*(D+1)/2 + (K-1); bic_values(i) = -2 * loglik + param_count * log(N); loglik_values(i) = loglik; end [~, best_idx] = min(bic_values); best_K = K_list(best_idx);

选BIC最小的K。注意如果你用的是fitgmdist，也可以直接调gmm.BIC，但自写代码的好处是能把loglik曲线和每个分量的有效点数一起打印出来，方便判断选中的K是否稳定。

我再分享一个调GMM的习惯：每次训练完，先打印loglik曲线形状、每个分量的点数Nk、以及det(Sigma(:,:,k))这三个值。如果loglik曲线单调上升、每个分量Nk都大于总样本数的5%、而且行列式没有接近0的值，这次训练基本是稳的。这套检查花不了十秒，但能省掉后面大量的返工。

有段时间我一直在折腾“分量数自适应”的方案，想着能不能让模型自动增减分量，后来发现BIC配合多次初始化已经足够解决绝大部分应用需求，自适应的复杂度反而引入更多不确定性。所以我现在更推荐把这个方向留给确有必要的人去研究，常规项目直接用BIC定K就够了。希望帮到你。

本文还有配套的精品资源，点击获取