拓冰网站定制 
简介这是一份GMM高斯混合聚类算法的Matlab完整实现基于matlab2023b编写适合计算机、电子信息工程、数学等专业学生用于课程设计、期末大作业和毕业设计也适合需要快速搭建聚类实验的研究者。代码采用参数化编程思路注释清晰结构紧凑便于二次开发聚类数目、协方差类型等关键参数集中可调运行后可直接输出聚类结果图以及不同协方差矩阵类型对模型性能的影响图省去繁琐调试一键即可复现典型实验。压缩包共2个文件包含1个主程序m脚本和1份特征数据xls表格整体仅17KB轻量易用。目前已有123人学习下载。借助这套源码读者可直观理解高斯混合模型在多簇数据上的建模过程、EM迭代与后验概率计算并通过协方差对比把握模型复杂度与拟合效果之间的权衡xls数据也能便捷替换为自有特征集用于课程实验或算法对比是一份教学与实战兼顾的参考实现。1. GMM聚类不是“高级K-means”而是一个概率生成模型拿到这份Matlab源码时我最先做的是直接运行main.m。特征数据.xls被读入后代码自动完成标准化、模型拟合和出图输出聚类结果图和协方差矩阵类型对模型性能的影响图。在matlab2023b环境下确实能做到一键出图这对课程设计、期末大作业来说非常省事但如果你只把它当成一个能输出标签的黑盒遇到重叠簇和不均衡数据时就会无从下手。GMM和K-means最大的区别在于GMM为每个样本计算出属于每个簇的后验概率而不是一个硬标签它假设数据由K个高斯成分加权混合生成。理解了这个概率生成视角再看后面的EM迭代、协方差矩阵类型和BIC选择才不会觉得参数设置是玄学。2. 高斯混合模型与EM算法先搞清楚GMM在算什么2.1 为什么单个高斯分布不够用假设特征只有两维数据自然地聚集在三个区域每个区域的散点大致呈椭圆形。用一个二元高斯分布去拟合得到的是一个中心点和一条协方差椭圆显然无法描述三个聚集中心。高斯混合模型把概率密度写成K个高斯成分的加权和p(x) sum_{k1}^{K} π_k · N(x | μ_k, Σ_k)其中π_k是第k个成分的混合权重满足π_k≥0且Σπ_k1。每个成分就是一个局部高斯分布它的均值μ_k决定了簇中心协方差矩阵Σ_k决定了该簇的形状和朝向。这种加权叠加让GMM能够拟合任意光滑的多峰分布这也是它常被用作密度估计器的原因。2.2 EM的E步和M步责任度与参数更新直接最大化对数似然没有闭合解因为每个样本具体来自哪个成分是隐藏的。EM算法通过两步交替迭代逼近最优参数E步固定当前参数计算每个样本x_i由第k个成分生成的后验概率γ_{ik}也叫责任度M步用这些责任度作为软权重重新估计π、μ、Σ。下面的Matlab风格代码展示了核心更新逻辑% X : N x d 数据矩阵K : 成分个数 % gamma: N x K 责任度矩阵 % E步计算每个样本对每个成分的责任度 log_pdf zeros(N, K); for k 1:K % mvnpdf 为多元高斯密度值 log_pdf(:, k) log(pi(k)) log(mvnpdf(X, mu(k,:), Sigma(:,:,k))); end % 减去每行最大值防止 exp 下溢 log_pdf log_pdf - max(log_pdf, [], 2); gamma exp(log_pdf) ./ sum(exp(log_pdf), 2); % M步按责任度加权更新参数 Nk sum(gamma, 1); % 每个成分的有效样本数 pi_new Nk / N; % 权重更新 for k 1:K mu_new(k,:) sum(gamma(:,k) .* X, 1) / Nk(k); X_centered X - mu_new(k,:); Sigma_new(:,:,k) (X_centered * (gamma(:,k) .* X_centered)) / Nk(k); end这段代码不是让你替换fitgmdist而是解释内置函数在迭代什么。E步里的数值稳定处理很重要直接计算exp(log_pdf)在数据维度高或者密度值很小时会下溢成0减掉每行最大值后结果不变但计算范围安全很多。M步本质上是加权平均责任度大的样本对μ和Σ的更新贡献更大。实际开发中我会用Matlab的fitgmdist直接拟合自己写EM还要额外处理协方差奇异、初始化不稳定等问题性价比太低。2.3 协方差矩阵类型与SharedCovariance参数的四种组合fitgmdist通过两个参数一起控制协方差结构CovarianceType决定每个成分的协方差是否只保留对角元素SharedCovariance决定所有成分是否共享同一个协方差矩阵。两两组合得到四种模型对应从“最灵活”到“最精简”的复杂度谱系CovarianceTypeSharedCovariance含义参数数量适用场景fullfalse每个成分独立的完整协方差矩阵最多各个簇形状、朝向差异大fulltrue所有成分共享同一个完整协方差矩阵中等各簇大小和朝向接近diagonalfalse每个成分独立但特征之间不相关中等偏少特征相关性弱要防过拟合diagonaltrue所有成分共享一个对角协方差矩阵最少样本少或只想粗略拟合这里要注意Matlab中没有“tied”这个选项共享协方差通过SharedCovariance为true实现别和Python的sklearn混淆。对二维特征来说full协方差能旋转椭圆方向而diagonal只能生成与坐标轴对齐的椭圆。原始代码里那张“协方差矩阵类型对模型性能的影响图”本质上就是在循环这四种组合后对比BIC或正确率。选择原则是样本量足够且簇形状差异明显时用fullfalse样本少或特征维度高时优先考虑diagonaltrue防止过拟合。2.4 用BIC/AIC给模型复杂度上锁GMM需要预先指定K但K通常未知。BIC和AIC是常见的选择依据BIC -2·logL p·log(N)AIC -2·logL 2·p其中p是参数个数。BIC对参数数量的惩罚更重适合样本量中等的情况。在Matlab中直接调用gmm.BIC可以拿到拟合模型的BIC值。如果想比较不同K可以这样写K_list 1:6; bic_values zeros(size(K_list)); gmm_models cell(size(K_list)); for i 1:length(K_list) gmm_models{i} fitgmdist(X, K_list(i), ... CovarianceType, full, ... SharedCovariance, false, ... Options, statset(Display, off)); bic_values(i) gmm_models{i}.BIC; end % 绘制 BIC 随 K 变化的曲线 figure; plot(K_list, bic_values, o-); xlabel(Number of components K); ylabel(BIC);这段循环把不同K下的模型存起来方便后续对比轮廓系数或直接导出标签。注意Fitting过程中如果某个K下出现协方差奇异fitgmdist会报错常见处理是给fitgmdist加RegularizationValue比如1e-6这在后面会提到。BIC越小代表在拟合度和复杂度之间取得平衡越好但不要只看BIC的最小值还要看曲线下降的“肘部”下降突然变缓的点往往比全局最小值更稳健。3. main.m源码逐段拆解从readmatrix到一键出图3.1 数据读取与标准化特征数据.xls进内存源码包里的main.m是整个流程的入口文件结构很清晰主要分为四个区段。第一段是数据读取与预处理。由于数据存放在特征数据.xls中源码使用readmatrix读取这一步比旧版xlsread更快且能自动识别数值列。我摘出的关键代码如下%% 数据读取与标准化 [data, ~, ~] xlsread(特征数据.xls); % 兼容旧版本也可用 readmatrix X data(:, 2:end); % 假设第一列是样本编号特征从第二列开始 % 如果数据没有编号列直接用 data 即可下面一行请按实际表头调整 % X data; % 标准化每个特征零均值、单位方差 [X_std, mu_X, sigma_X] zscore(X);readmatrix在R2023b中是官方推荐的数据导入函数能处理.xls和.xlsx不依赖Excel插件。这里假设第一列是编号如果原始文件没有这一列就去掉这行改成X data;。zscore按列做标准化把不同量纲的特征拉到同一尺度同时返回每个原始特征的均值mu_X和标准差sigma_X供后面新样本预测时复用。在GMM里特征尺度会直接影响协方差矩阵的奇异性比如一个特征范围在0~1另一个在0~10000未经标准化直接拟合协方差矩阵很可能接近奇异导致EM迭代失败。3.2 fitgmdist核心参数与聚类输出标准化之后进入模型拟合源码中通过参数化编程把K、协方差类型、展示开关都集中定义在文件头部方便修改。核心调用如下%% GMM 拟合 K 3; % 预设聚类数可按业务调整 covType full; % full 或 diagonal shareCov false; % true 或 false regVal 1e-6; % 正则化项防止协方差奇异 gmm fitgmdist(X_std, K, ... CovarianceType, covType, ... SharedCovariance, shareCov, ... RegularizationValue, regVal, ... Options, statset(Display, final)); % 硬分类标签与软分类概率 [label, ~, posterior] cluster(gmm, X_std);fitgmdist的参数需要逐一说明参数名可选值作用K正整数混合成分数量即聚类数CovarianceTypefull / diagonal是否只保留对角线协方差SharedCovariancetrue / false是否所有成分共享协方差RegularizationValue如1e-6防止协方差奇异Optionsstatset(...)控制迭代显示和容差cluster函数返回三个输出这里用~跳过第二个负对数似然nlogl保留label和posterior。label是每个样本最大后验概率对应的类别编号posterior是N×K的概率矩阵其中每一行和为1。后续绘制聚类图我们用的是label而异常检测则直接利用posterior。3.3 聚类结果图的绘制逻辑一份能直接放进课程设计报告的图至少要包含聚类散点和聚类中心。源码用gscatter按label绘制不同颜色的散点再叠加GMM的均值位置示例如下%% 聚类结果图 figure(Color, w); gscatter(X_std(:,1), X_std(:,2), label, rgbkcmy, ., 15); hold on; plot(gmm.mu(:,1), gmm.mu(:,2), kx, MarkerSize, 12, LineWidth, 2); xlabel(Feature 1 (Standardized)); ylabel(Feature 2 (Standardized)); title(GMM Clustering Result (K3)); legend(Cluster 1, Cluster 2, Cluster 3, Centroid); grid on;gscatter是Matlab的Grouped Scatter要求X_std至少有两列高维特征时通常会取PCA前两维投影后再画。我一般会先对高维数据做pca再用前两个主成分替换X_std输入gscatter这样不会失真太严重。plot里的gmm.mu是K×d矩阵直接取前两列标出各成分中心。如果特征维度不是2这张图只能反映二维投影报告里需要注明是前两维或主成分。3.4 协方差矩阵影响对比图的批量生成源码的第二张输出图是“协方差矩阵类型对模型性能的影响图”它循环四种协方差组合拟合并记录BIC和迭代时间最后用bar或plot对比。核心结构如下%% 对比四种协方差设置 covSettings {full, full, diagonal, diagonal}; shareSettings [false, true, false, true]; bicList zeros(1, 4); for i 1:4 gmm_temp fitgmdist(X_std, K, ... CovarianceType, covSettings{i}, ... SharedCovariance, shareSettings(i), ... RegularizationValue, 1e-6, ... Options, statset(Display, off)); bicList(i) gmm_temp.BIC; end figure(Color, w); bar(categorical({full独立,full共享,diag独立,diag共享}), bicList); ylabel(BIC); title(Effect of Covariance Type on GMM);这里categorical把四个组合的名称作为横轴标签让BIC对比一目了然。从工程角度看这张图的价值在于如果BIC差异很小优先选参数少的diagonalsharedTrue因为它泛化更好训练也更快。如果BIC差异巨大说明协方差结构对数据集敏感需要结合业务判断是哪些簇的形状差异导致的。运行这段代码时如果某个组合不收敛别急着加正则先检查X_std是否标准化以及K是否设得过大。4. 协方差类型、初始化与聚类数K四个影响聚类效果的坑4.1 四种协方差组合的适用场景上一章的对比图能展示哪种协方差结构在数据上BIC最低但实际应用中还需要结合样本量和业务语义。下面这张表是我在使用这类Matlab源码时的参考数据特征推荐组合原因样本量大各簇大小形状不同full sharedfalse最能表示簇的真实边界样本量中等簇接近球形或椭圆方向一致full sharedtrue降低方差避免过拟合特征之间有较强相关性但样本量少diagonal sharedfalse每簇独立却不需要估计协方差中非对角项样本量极少特征独立diagonal sharedtrue参数最少最稳定这个表格直接对应main.m中的covType和shareCov变量。如果你的实验结果是diagsharedtrue的BIC最低不要觉得意外那说明你的特征之间相关性弱或者数据量不足以支撑完整协方差矩阵的估计。相反如果fullfalse的BIC显著更优说明簇的形状差异真实存在这时候强行用diagonal会得到边缘模糊的聚类边界。4.2 初始化与Replicates的坑fitgmdist默认使用k-means的变体初始化但EM算法对初始值敏感容易陷入局部最优。即使定了K和协方差结构一次运行的结果也可能不同。解决办法是设置Replicates参数让Matlab从多个随机起点拟合选出最大似然的模型gmm fitgmdist(X_std, 3, ... CovarianceType, full, ... SharedCovariance, false, ... Replicates, 10, ... % 10次随机初始化保留最优 Options, statset(MaxIter, 500));Replicates每多一次运行时间大约翻一倍但能显著降低局部最优的风险。注意当运行fitgmdist报错提示协方差矩阵不正定时优先检查K和RegularizationValue而不是盲目加大Replicates。另一个更实际的坑是如果K过大某个成分可能只有一两个样本它的协方差矩阵会退化。这时RegularizationValue要适当加大比如从1e-6提到1e-4否则fitgmdist直接报错。遇到这种错误先检查K是否明显大于类簇数再考虑加正则。4.3 聚类数K的确定BIC、轮廓系数与稳定性上一章对比图只针对固定K但K的选择才是GMM最关键的决策。我用一个循环同时计算BIC和轮廓系数综合判断K_list 1:6; silList zeros(size(K_list)); for i 1:length(K_list) gmm_temp fitgmdist(X_std, K_list(i), ... CovarianceType, full, ... SharedCovariance, false, ... RegularizationValue, 1e-6, ... Replicates, 5); [tmp_label, ~] cluster(gmm_temp, X_std); if K_list(i) 1 silList(i) 0; % 单簇没有轮廓系数 else silList(i) mean(silhouette(X_std, tmp_label)); end end轮廓系数衡量簇内紧密度和簇间分离度Matlab的silhouette函数返回每个点的轮廓值取平均得到整体指标。BIC最小值并不总是对应业务上可解释的划分轮廓系数倾向于紧致且分离的簇两者结合能避免被单一指标带偏。如果K3时BIC曲线已明显变平而K5时轮廓系数最大我一般会选择K3因为简单模型在工程上更稳定。4.4 先标准化还是先聚类特征尺度与异常值很多直接跑源码的初学者会跳过zscore结果聚类图全被某个量纲大的特征主导。GMM对尺度敏感因为它依赖欧氏距离计算协方差矩阵。zscore将每个特征变成均值0方差1是一种公平处理但也抹掉了特征间的相对重要性。如果业务上知道某特征更重要可以用自定义权重替代标准化比如w [1, 3, 1]; % 给第二个特征三倍权重 X_weighted X .* w; % 乘性加权等价于调整协方差贡献这里更推荐的加权方式是对每个特征除以设定的尺度参数而不是简单乘权重。异常值也是坑单个离群点会把某个成分的均值拉偏协方差矩阵被放大聚类结果整体漂移。我会先画箱线图或使用robustcov检查如果异常值超过1%先做截尾处理再进GMM。对于课程设计这些预处理步骤最好单独放在main.m文件头部并加注释保证评估老师能一眼看到。5. 把GMM当密度估计器异常检测与软分类的进阶用法5.1 用后验概率做软分类聚类结果图中label给出的是硬标签但GMM真正有价值的输出是后验概率。比如新来了一个样本x_new可以用posterior函数计算它属于每个成分的概率然后取最大者作为软分类结果。这在用户分群场景中很有用可能某个用户同时有70%概率属于高活跃群、30%属于中等活跃群硬分类会丢掉这层信息。代码上就是x_new [0.5, -1.2]; % 新样本需与训练时同样的标准化 x_new_std (x_new - mu_X) ./ sigma_X; % 用训练集均值和标准差标准化 [poster_new, ~] posterior(gmm, x_new_std); [prob, idx] max(poster_new); fprintf(样本属于第%d个簇概率为%.2f\n, idx, prob);注意训练时的mu_X、sigma_X是在zscore前用mean和std计算出来的在3.1节已经通过zscore的第二个、第三个输出保存。后验概率的可视化可以画成热度图如果特征大于三维就画两两投影矩阵这篇文章里不再展开。5.2 基于对数似然的异常检测阈值当GMM拟合完成后它可以作为一个密度估计器来用。每个样本的对数似然log(p(x))越低说明它越不可能来自这个生成过程越可能是异常点。main.m里虽然没直接提供异常检测代码但只需一行调用即可loglike log(pdf(gmm, X_std)); % N x 1 的对数密度值 threshold prctile(loglike, 5); % 取5%分位点作为异常阈值 outlier_idx find(loglike threshold);这个阈值取5%是比较常见的默认做法但具体分位点要看业务。如果数据集本身干净用低于最小值的5%会误杀正常点如果数据确实包含大量噪点可以把阈值上移到10%。这里有个验证技巧把异常点标红画到聚类图上如果它们都落在椭圆边界外说明GMM的密度估计和聚类几何是一致的如果异常点散布在簇内部则说明数据分布不适合用单一高斯成分建模需要考虑增加K或改用其他密度估计模型。本文还有配套的精品资源点击获取 GMM高斯混合模型聚类原理与Matlab实现：从EM算法到BIC选择
简介这是一份GMM高斯混合聚类算法的Matlab完整实现基于matlab2023b编写适合计算机、电子信息工程、数学等专业学生用于课程设计、期末大作业和毕业设计也适合需要快速搭建聚类实验的研究者。代码采用参数化编程思路注释清晰结构紧凑便于二次开发聚类数目、协方差类型等关键参数集中可调运行后可直接输出聚类结果图以及不同协方差矩阵类型对模型性能的影响图省去繁琐调试一键即可复现典型实验。压缩包共2个文件包含1个主程序m脚本和1份特征数据xls表格整体仅17KB轻量易用。目前已有123人学习下载。借助这套源码读者可直观理解高斯混合模型在多簇数据上的建模过程、EM迭代与后验概率计算并通过协方差对比把握模型复杂度与拟合效果之间的权衡xls数据也能便捷替换为自有特征集用于课程实验或算法对比是一份教学与实战兼顾的参考实现。1. GMM聚类不是“高级K-means”而是一个概率生成模型拿到这份Matlab源码时我最先做的是直接运行main.m。特征数据.xls被读入后代码自动完成标准化、模型拟合和出图输出聚类结果图和协方差矩阵类型对模型性能的影响图。在matlab2023b环境下确实能做到一键出图这对课程设计、期末大作业来说非常省事但如果你只把它当成一个能输出标签的黑盒遇到重叠簇和不均衡数据时就会无从下手。GMM和K-means最大的区别在于GMM为每个样本计算出属于每个簇的后验概率而不是一个硬标签它假设数据由K个高斯成分加权混合生成。理解了这个概率生成视角再看后面的EM迭代、协方差矩阵类型和BIC选择才不会觉得参数设置是玄学。2. 高斯混合模型与EM算法先搞清楚GMM在算什么2.1 为什么单个高斯分布不够用假设特征只有两维数据自然地聚集在三个区域每个区域的散点大致呈椭圆形。用一个二元高斯分布去拟合得到的是一个中心点和一条协方差椭圆显然无法描述三个聚集中心。高斯混合模型把概率密度写成K个高斯成分的加权和p(x) sum_{k1}^{K} π_k · N(x | μ_k, Σ_k)其中π_k是第k个成分的混合权重满足π_k≥0且Σπ_k1。每个成分就是一个局部高斯分布它的均值μ_k决定了簇中心协方差矩阵Σ_k决定了该簇的形状和朝向。这种加权叠加让GMM能够拟合任意光滑的多峰分布这也是它常被用作密度估计器的原因。2.2 EM的E步和M步责任度与参数更新直接最大化对数似然没有闭合解因为每个样本具体来自哪个成分是隐藏的。EM算法通过两步交替迭代逼近最优参数E步固定当前参数计算每个样本x_i由第k个成分生成的后验概率γ_{ik}也叫责任度M步用这些责任度作为软权重重新估计π、μ、Σ。下面的Matlab风格代码展示了核心更新逻辑% X : N x d 数据矩阵K : 成分个数 % gamma: N x K 责任度矩阵 % E步计算每个样本对每个成分的责任度 log_pdf zeros(N, K); for k 1:K % mvnpdf 为多元高斯密度值 log_pdf(:, k) log(pi(k)) log(mvnpdf(X, mu(k,:), Sigma(:,:,k))); end % 减去每行最大值防止 exp 下溢 log_pdf log_pdf - max(log_pdf, [], 2); gamma exp(log_pdf) ./ sum(exp(log_pdf), 2); % M步按责任度加权更新参数 Nk sum(gamma, 1); % 每个成分的有效样本数 pi_new Nk / N; % 权重更新 for k 1:K mu_new(k,:) sum(gamma(:,k) .* X, 1) / Nk(k); X_centered X - mu_new(k,:); Sigma_new(:,:,k) (X_centered * (gamma(:,k) .* X_centered)) / Nk(k); end这段代码不是让你替换fitgmdist而是解释内置函数在迭代什么。E步里的数值稳定处理很重要直接计算exp(log_pdf)在数据维度高或者密度值很小时会下溢成0减掉每行最大值后结果不变但计算范围安全很多。M步本质上是加权平均责任度大的样本对μ和Σ的更新贡献更大。实际开发中我会用Matlab的fitgmdist直接拟合自己写EM还要额外处理协方差奇异、初始化不稳定等问题性价比太低。2.3 协方差矩阵类型与SharedCovariance参数的四种组合fitgmdist通过两个参数一起控制协方差结构CovarianceType决定每个成分的协方差是否只保留对角元素SharedCovariance决定所有成分是否共享同一个协方差矩阵。两两组合得到四种模型对应从“最灵活”到“最精简”的复杂度谱系CovarianceTypeSharedCovariance含义参数数量适用场景fullfalse每个成分独立的完整协方差矩阵最多各个簇形状、朝向差异大fulltrue所有成分共享同一个完整协方差矩阵中等各簇大小和朝向接近diagonalfalse每个成分独立但特征之间不相关中等偏少特征相关性弱要防过拟合diagonaltrue所有成分共享一个对角协方差矩阵最少样本少或只想粗略拟合这里要注意Matlab中没有“tied”这个选项共享协方差通过SharedCovariance为true实现别和Python的sklearn混淆。对二维特征来说full协方差能旋转椭圆方向而diagonal只能生成与坐标轴对齐的椭圆。原始代码里那张“协方差矩阵类型对模型性能的影响图”本质上就是在循环这四种组合后对比BIC或正确率。选择原则是样本量足够且簇形状差异明显时用fullfalse样本少或特征维度高时优先考虑diagonaltrue防止过拟合。2.4 用BIC/AIC给模型复杂度上锁GMM需要预先指定K但K通常未知。BIC和AIC是常见的选择依据BIC -2·logL p·log(N)AIC -2·logL 2·p其中p是参数个数。BIC对参数数量的惩罚更重适合样本量中等的情况。在Matlab中直接调用gmm.BIC可以拿到拟合模型的BIC值。如果想比较不同K可以这样写K_list 1:6; bic_values zeros(size(K_list)); gmm_models cell(size(K_list)); for i 1:length(K_list) gmm_models{i} fitgmdist(X, K_list(i), ... CovarianceType, full, ... SharedCovariance, false, ... Options, statset(Display, off)); bic_values(i) gmm_models{i}.BIC; end % 绘制 BIC 随 K 变化的曲线 figure; plot(K_list, bic_values, o-); xlabel(Number of components K); ylabel(BIC);这段循环把不同K下的模型存起来方便后续对比轮廓系数或直接导出标签。注意Fitting过程中如果某个K下出现协方差奇异fitgmdist会报错常见处理是给fitgmdist加RegularizationValue比如1e-6这在后面会提到。BIC越小代表在拟合度和复杂度之间取得平衡越好但不要只看BIC的最小值还要看曲线下降的“肘部”下降突然变缓的点往往比全局最小值更稳健。3. main.m源码逐段拆解从readmatrix到一键出图3.1 数据读取与标准化特征数据.xls进内存源码包里的main.m是整个流程的入口文件结构很清晰主要分为四个区段。第一段是数据读取与预处理。由于数据存放在特征数据.xls中源码使用readmatrix读取这一步比旧版xlsread更快且能自动识别数值列。我摘出的关键代码如下%% 数据读取与标准化 [data, ~, ~] xlsread(特征数据.xls); % 兼容旧版本也可用 readmatrix X data(:, 2:end); % 假设第一列是样本编号特征从第二列开始 % 如果数据没有编号列直接用 data 即可下面一行请按实际表头调整 % X data; % 标准化每个特征零均值、单位方差 [X_std, mu_X, sigma_X] zscore(X);readmatrix在R2023b中是官方推荐的数据导入函数能处理.xls和.xlsx不依赖Excel插件。这里假设第一列是编号如果原始文件没有这一列就去掉这行改成X data;。zscore按列做标准化把不同量纲的特征拉到同一尺度同时返回每个原始特征的均值mu_X和标准差sigma_X供后面新样本预测时复用。在GMM里特征尺度会直接影响协方差矩阵的奇异性比如一个特征范围在0~1另一个在0~10000未经标准化直接拟合协方差矩阵很可能接近奇异导致EM迭代失败。3.2 fitgmdist核心参数与聚类输出标准化之后进入模型拟合源码中通过参数化编程把K、协方差类型、展示开关都集中定义在文件头部方便修改。核心调用如下%% GMM 拟合 K 3; % 预设聚类数可按业务调整 covType full; % full 或 diagonal shareCov false; % true 或 false regVal 1e-6; % 正则化项防止协方差奇异 gmm fitgmdist(X_std, K, ... CovarianceType, covType, ... SharedCovariance, shareCov, ... RegularizationValue, regVal, ... Options, statset(Display, final)); % 硬分类标签与软分类概率 [label, ~, posterior] cluster(gmm, X_std);fitgmdist的参数需要逐一说明参数名可选值作用K正整数混合成分数量即聚类数CovarianceTypefull / diagonal是否只保留对角线协方差SharedCovariancetrue / false是否所有成分共享协方差RegularizationValue如1e-6防止协方差奇异Optionsstatset(...)控制迭代显示和容差cluster函数返回三个输出这里用~跳过第二个负对数似然nlogl保留label和posterior。label是每个样本最大后验概率对应的类别编号posterior是N×K的概率矩阵其中每一行和为1。后续绘制聚类图我们用的是label而异常检测则直接利用posterior。3.3 聚类结果图的绘制逻辑一份能直接放进课程设计报告的图至少要包含聚类散点和聚类中心。源码用gscatter按label绘制不同颜色的散点再叠加GMM的均值位置示例如下%% 聚类结果图 figure(Color, w); gscatter(X_std(:,1), X_std(:,2), label, rgbkcmy, ., 15); hold on; plot(gmm.mu(:,1), gmm.mu(:,2), kx, MarkerSize, 12, LineWidth, 2); xlabel(Feature 1 (Standardized)); ylabel(Feature 2 (Standardized)); title(GMM Clustering Result (K3)); legend(Cluster 1, Cluster 2, Cluster 3, Centroid); grid on;gscatter是Matlab的Grouped Scatter要求X_std至少有两列高维特征时通常会取PCA前两维投影后再画。我一般会先对高维数据做pca再用前两个主成分替换X_std输入gscatter这样不会失真太严重。plot里的gmm.mu是K×d矩阵直接取前两列标出各成分中心。如果特征维度不是2这张图只能反映二维投影报告里需要注明是前两维或主成分。3.4 协方差矩阵影响对比图的批量生成源码的第二张输出图是“协方差矩阵类型对模型性能的影响图”它循环四种协方差组合拟合并记录BIC和迭代时间最后用bar或plot对比。核心结构如下%% 对比四种协方差设置 covSettings {full, full, diagonal, diagonal}; shareSettings [false, true, false, true]; bicList zeros(1, 4); for i 1:4 gmm_temp fitgmdist(X_std, K, ... CovarianceType, covSettings{i}, ... SharedCovariance, shareSettings(i), ... RegularizationValue, 1e-6, ... Options, statset(Display, off)); bicList(i) gmm_temp.BIC; end figure(Color, w); bar(categorical({full独立,full共享,diag独立,diag共享}), bicList); ylabel(BIC); title(Effect of Covariance Type on GMM);这里categorical把四个组合的名称作为横轴标签让BIC对比一目了然。从工程角度看这张图的价值在于如果BIC差异很小优先选参数少的diagonalsharedTrue因为它泛化更好训练也更快。如果BIC差异巨大说明协方差结构对数据集敏感需要结合业务判断是哪些簇的形状差异导致的。运行这段代码时如果某个组合不收敛别急着加正则先检查X_std是否标准化以及K是否设得过大。4. 协方差类型、初始化与聚类数K四个影响聚类效果的坑4.1 四种协方差组合的适用场景上一章的对比图能展示哪种协方差结构在数据上BIC最低但实际应用中还需要结合样本量和业务语义。下面这张表是我在使用这类Matlab源码时的参考数据特征推荐组合原因样本量大各簇大小形状不同full sharedfalse最能表示簇的真实边界样本量中等簇接近球形或椭圆方向一致full sharedtrue降低方差避免过拟合特征之间有较强相关性但样本量少diagonal sharedfalse每簇独立却不需要估计协方差中非对角项样本量极少特征独立diagonal sharedtrue参数最少最稳定这个表格直接对应main.m中的covType和shareCov变量。如果你的实验结果是diagsharedtrue的BIC最低不要觉得意外那说明你的特征之间相关性弱或者数据量不足以支撑完整协方差矩阵的估计。相反如果fullfalse的BIC显著更优说明簇的形状差异真实存在这时候强行用diagonal会得到边缘模糊的聚类边界。4.2 初始化与Replicates的坑fitgmdist默认使用k-means的变体初始化但EM算法对初始值敏感容易陷入局部最优。即使定了K和协方差结构一次运行的结果也可能不同。解决办法是设置Replicates参数让Matlab从多个随机起点拟合选出最大似然的模型gmm fitgmdist(X_std, 3, ... CovarianceType, full, ... SharedCovariance, false, ... Replicates, 10, ... % 10次随机初始化保留最优 Options, statset(MaxIter, 500));Replicates每多一次运行时间大约翻一倍但能显著降低局部最优的风险。注意当运行fitgmdist报错提示协方差矩阵不正定时优先检查K和RegularizationValue而不是盲目加大Replicates。另一个更实际的坑是如果K过大某个成分可能只有一两个样本它的协方差矩阵会退化。这时RegularizationValue要适当加大比如从1e-6提到1e-4否则fitgmdist直接报错。遇到这种错误先检查K是否明显大于类簇数再考虑加正则。4.3 聚类数K的确定BIC、轮廓系数与稳定性上一章对比图只针对固定K但K的选择才是GMM最关键的决策。我用一个循环同时计算BIC和轮廓系数综合判断K_list 1:6; silList zeros(size(K_list)); for i 1:length(K_list) gmm_temp fitgmdist(X_std, K_list(i), ... CovarianceType, full, ... SharedCovariance, false, ... RegularizationValue, 1e-6, ... Replicates, 5); [tmp_label, ~] cluster(gmm_temp, X_std); if K_list(i) 1 silList(i) 0; % 单簇没有轮廓系数 else silList(i) mean(silhouette(X_std, tmp_label)); end end轮廓系数衡量簇内紧密度和簇间分离度Matlab的silhouette函数返回每个点的轮廓值取平均得到整体指标。BIC最小值并不总是对应业务上可解释的划分轮廓系数倾向于紧致且分离的簇两者结合能避免被单一指标带偏。如果K3时BIC曲线已明显变平而K5时轮廓系数最大我一般会选择K3因为简单模型在工程上更稳定。4.4 先标准化还是先聚类特征尺度与异常值很多直接跑源码的初学者会跳过zscore结果聚类图全被某个量纲大的特征主导。GMM对尺度敏感因为它依赖欧氏距离计算协方差矩阵。zscore将每个特征变成均值0方差1是一种公平处理但也抹掉了特征间的相对重要性。如果业务上知道某特征更重要可以用自定义权重替代标准化比如w [1, 3, 1]; % 给第二个特征三倍权重 X_weighted X .* w; % 乘性加权等价于调整协方差贡献这里更推荐的加权方式是对每个特征除以设定的尺度参数而不是简单乘权重。异常值也是坑单个离群点会把某个成分的均值拉偏协方差矩阵被放大聚类结果整体漂移。我会先画箱线图或使用robustcov检查如果异常值超过1%先做截尾处理再进GMM。对于课程设计这些预处理步骤最好单独放在main.m文件头部并加注释保证评估老师能一眼看到。5. 把GMM当密度估计器异常检测与软分类的进阶用法5.1 用后验概率做软分类聚类结果图中label给出的是硬标签但GMM真正有价值的输出是后验概率。比如新来了一个样本x_new可以用posterior函数计算它属于每个成分的概率然后取最大者作为软分类结果。这在用户分群场景中很有用可能某个用户同时有70%概率属于高活跃群、30%属于中等活跃群硬分类会丢掉这层信息。代码上就是x_new [0.5, -1.2]; % 新样本需与训练时同样的标准化 x_new_std (x_new - mu_X) ./ sigma_X; % 用训练集均值和标准差标准化 [poster_new, ~] posterior(gmm, x_new_std); [prob, idx] max(poster_new); fprintf(样本属于第%d个簇概率为%.2f\n, idx, prob);注意训练时的mu_X、sigma_X是在zscore前用mean和std计算出来的在3.1节已经通过zscore的第二个、第三个输出保存。后验概率的可视化可以画成热度图如果特征大于三维就画两两投影矩阵这篇文章里不再展开。5.2 基于对数似然的异常检测阈值当GMM拟合完成后它可以作为一个密度估计器来用。每个样本的对数似然log(p(x))越低说明它越不可能来自这个生成过程越可能是异常点。main.m里虽然没直接提供异常检测代码但只需一行调用即可loglike log(pdf(gmm, X_std)); % N x 1 的对数密度值 threshold prctile(loglike, 5); % 取5%分位点作为异常阈值 outlier_idx find(loglike threshold);这个阈值取5%是比较常见的默认做法但具体分位点要看业务。如果数据集本身干净用低于最小值的5%会误杀正常点如果数据确实包含大量噪点可以把阈值上移到10%。这里有个验证技巧把异常点标红画到聚类图上如果它们都落在椭圆边界外说明GMM的密度估计和聚类几何是一致的如果异常点散布在簇内部则说明数据分布不适合用单一高斯成分建模需要考虑增加K或改用其他密度估计模型。本文还有配套的精品资源点击获取