简介高斯混合模型GMM与高斯混合回归GMR是机器学习中常用的概率建模工具前者用于无监督聚类和密度估计后者在此基础上实现连续变量的回归预测。本资源提供一套完整的MATLAB实现代码目标使用者是具备基础统计知识、正在学习机器学习或需要快速搭建GMM/GMR模型的研究人员和学生。压缩包内共26个文件其中19个.m脚本覆盖了高斯概率计算、EM算法中的E步与M步、K均值初始化、混合回归训练与预测、LQR控制等核心流程6个.mat数据文件用于直接运行和复现演示另有1个txt说明文档便于理清代码结构。资源包大小仅98KB轻量便携。目前已有1958人学习代码中配有实例和可视化图形通过运行demo可直观理解参数估计迭代过程和回归效果。借助该资源可以系统掌握GMM聚类的数学原理、GMR回归的预测逻辑以及MATLAB工程化实现思路是兼顾理论学习与实际建模的优质参考。 做数据分析的同行应该都有过这种经历数据分布明明有好几个山头硬拿一个高斯分布去拟合估计出来的均值落在两座山之间的谷底方差也大得离谱完全不能代表真实结构。这时候就得搬出高斯混合模型GMM。用K个高斯分布加权叠加去逼近任意复杂的多峰密度聚类、异常检测、密度估计、回归预测量它都能干。而高斯混合回归则是把GMM从描述数据升级成预测输出在软切换的局部线性回归任务里比单纯套一条全局直线要实用得多。这篇文章不绕弯子直接从原理讲到MATLAB实现把GMM、EM迭代、高斯混合回归的预测函数、常见坑位一次说清楚适合正在做模式识别、信号处理、系统辨识和机器学习相关工作的同学参考。1. GMM到底在算什么从单峰到多峰的概率视角1.1 一个高斯分布为什么不够用高斯分布是统计建模里最常用的假设均值和协方差两个东西就能完整描述一个对称的钟形分布。但实际问题里数据往往不是一个整体。比如同一个工艺线上生产出来的零件可能来自两台精度不同的机床传感器采集的振动信号可能包含正常工况和异常工况两种状态。把这些数据混在一起看直方图大概率会出现两个甚至多个峰。如果用单一高斯分布去拟合这种多峰数据最大似然估计出来的均值会把所有峰值揉在一起结果是既不在这也不在那。协方差也会被迫拉大覆盖两个峰之间的空白区域。密度估计严重失真后续的检测、分类都会跟着出错。这时候把数据拆成几个高斯分量的叠加让每个分量负责描述一个局部聚集结构再按权重组合起来就能得到非常接近真实分布的密度函数——这就是GMM存在的意义。1.2 高斯混合模型的数学表达GMM的密度函数形式上很干净[ p(x)\sum_{k1}^{K}\pi_k,\mathcal{N}(x|\mu_k,\Sigma_k) ]其中 (K) 是分量个数(\pi_k) 是第 (k) 个分量的混合系数满足 (\pi_k\ge 0) 且 (\sum_k\pi_k1)(\mu_k) 和 (\Sigma_k) 分别是第 (k) 个高斯分量的均值向量和协方差矩阵。可以把 (\pi_k) 理解为数据里大约有 (\pi_k) 比例的点是由第 (k) 个分量产生的。举个最简单的例子两个分量(\pi_10.6)(\pi_20.4)第一条高斯中心在 ((1,1))第二条中心在 ((5,5))那么最终生成的数据会呈现两个椭圆簇一个密一些一个疏一些。这里的隐藏变量就是每个样本到底属于哪个分量这个变量我们观测不到所以才需要EM算法来迭代求解。1.3 什么时候优先选GMMGMM不是万能的但它有几个特别合适的场景密度估计面对多峰分布需要给任意一个点输出概率密度值。软聚类像K-means那种硬分配每个点只能属于一个簇GMM给出的是属于每个簇的后验概率对重叠簇更友好。异常检测用训练好的GMM计算样本的似然概率概率低于阈值就判为异常。作为回归或降维的前置模块高斯混合回归、高斯混合隐变量模型等都是在GMM的基本框架上加任务。如果只是做大规模硬聚类K-means更快更省事但当你需要概率、需要椭圆形的簇边界、需要把聚类结果进一步接回归或分类时GMM是更合适的底座。2. EM算法GMM参数估计的核心引擎2.1 为什么极大似然直接做不动GMM的参数包括 (\pi_k,\mu_k,\Sigma_k)。如果直接用极大似然估计需要对数似然函数[ \log p(X)\sum_{i1}^{N}\log\left(\sum_{k1}^{K}\pi_k,\mathcal{N}(x_i|\mu_k,\Sigma_k)\right) ]问题就出在 (\log) 里面套了一个求和。求导之后混合系数的更新会和所有未知参数纠缠在一起得不到闭式解。于是EM算法出场先假设参数已知估计每个样本的隐变量再把隐变量当作已知重新估计参数。两步交替不断抬升对数似然。2.2 E步算责任值在第 (t) 次迭代中用当前参数 ({\pi_k^{(t)},\mu_k^{(t)},\Sigma_k^{(t)}}) 计算第 (i) 个样本属于第 (k) 个分量的后验概率也叫责任值[ \gamma_{ik}\frac{\pi_k^{(t)},\mathcal{N}(x_i|\mu_k^{(t)},\Sigma_k^{(t)})}{\sum_{j1}^{K}\pi_j^{(t)},\mathcal{N}(x_i|\mu_j^{(t)},\Sigma_j^{(t)})} ]这个式子的分母是所有分量的加权密度之和归一化后(\gamma_{ik}) 表示样本 (x_i) 有多大概率来自分量 (k)。如果在某个分量附近责任值就趋近1如果落在两个分量中间责任值会按两边的密度占比分配。2.3 M步按责任值加权更新参数拿到责任值后参数更新变得很直观[ \pi_k^{(t1)}\frac{1}{N}\sum_{i1}^{N}\gamma_{ik} ][ \mu_k^{(t1)}\frac{\sum_{i1}^{N}\gamma_{ik}x_i}{\sum_{i1}^{N}\gamma_{ik}} ][ \Sigma_k^{(t1)}\frac{\sum_{i1}^{N}\gamma_{ik}(x_i-\mu_k^{(t1)})(x_i-\mu_k^{(t1)})^T}{\sum_{i1}^{N}\gamma_{ik}} ]本质上就是加权平均。责任值大的样本对当前分量的均值和协方差贡献大责任值小的贡献小。迭代交替进行直到对数似然的增量小于某个阈值或者达到最大迭代次数。我手写EM时初始化一般用K-means的聚类结果把K-means得到的簇心作为 (\mu) 的初值按簇内样本协方差作为 (\Sigma) 初值(\pi) 设为各类样本占比。这样比随机初始化收敛稳定得多也能有效避开一部分局部最优。3. MATLAB实现GMM现成函数与手写EM两条路3.1 最省事的方式fitgmdistMATLAB里做GMM第一选择是统计工具箱自带的fitgmdist。它内部已经把初始化、EM迭代、协方差类型选择、正则化都封装好了。基本用法rng(42); % 生成两个簇的数据 X [randn(300, 2) * 0.6 [1 1]; randn(200, 2) * 1.2 [5 5]]; k 2; gm fitgmdist(X, k, ... CovarianceType, full, ... RegularizationValue, 1e-6);fitgmdist返回的gm是gmdistribution对象直接用pdf(gm, X)算密度用cluster(gm, X)做硬聚类用random(gm, n)生成服从该混合模型的样本。CovarianceType可以选full或diagonal前者能拟合椭圆形的簇后者计算更稳更快RegularizationValue会在协方差矩阵对角线上加一个小常数防止矩阵奇异这个参数在样本少、维度高时特别关键。3.2 手写一个简化版EM只调现成函数容易踩到不知道内部发生了什么的坑。下面这段代码是简化版EM适合学习function [mu, Sigma, pi_, llh] my_gmm_em(X, K, maxIter, tol) % 输入X为N行D列的数据K为分量数 % 输出mu为K行D列Sigma为D×D×K数组pi_为K维向量 N size(X, 1); D size(X, 2); % 用kmeans做初始化 [idx, C] kmeans(X, K, Replicates, 5); mu C; pi_ zeros(K, 1); Sigma zeros(D, D, K); for k 1:K Xk X(idx k, :); pi_(k) size(Xk, 1) / N; Sigma(:, :, k) cov(Xk) 1e-6 * eye(D); end llh zeros(maxIter, 1); for iter 1:maxIter % E步 resp zeros(N, K); for k 1:K resp(:, k) pi_(k) * mvnpdf(X, mu(k, :), Sigma(:, :, k)); end sum_resp sum(resp, 2); llh(iter) sum(log(sum_resp)); resp resp ./ sum_resp; % 归一化成责任值 % 检查收敛 if iter 1 abs(llh(iter) - llh(iter-1)) tol llh llh(1:iter); break; end % M步 Nk sum(resp, 1); % 1×K pi_ Nk / N; mu (resp * X) ./ Nk; % K×D for k 1:K Xc X - mu(k, :); Sigma(:, :, k) (Xc .* resp(:, k)) * Xc / Nk(k) 1e-6 * eye(D); end end end这个版本的M步里Xc .* resp(:, k)用到了MATLAB的隐式扩展Xc是 (N\times D) 的零均值样本矩阵resp(:, k)是 (N\times 1) 的权重列向量乘出来之后每一行都乘上对应样本的责任值。再左乘Xc就是在计算加权外积和最后除以 (N_k)。如果你用旧版MATLAB可能不支持隐式扩展可以改成repmat。3.3 聚类结果可视化硬聚类只用一行idx cluster(gm, X); gscatter(X(:, 1), X(:, 2), idx);如果想画每个分量的置信椭圆需要从gm.mu和gm.Sigma里取均值与协方差再在椭圆方程里算角度和半轴长度。一个常用思路是对协方差矩阵做特征值分解特征值开方后乘以对应方向的单位向量就得到椭圆的长短轴。配合gscatter出来的图会专业很多。4. 高斯混合回归从联合分布到条件预测4.1 核心思路把y也当成高斯分量的一部分高斯混合回归Gaussian Mixture Regression简称GMR和普通回归的建模路径不一样。普通回归直接建模 (p(y|x))而GMR先把输入 (x) 和输出 (y) 拼成一个向量对联合分布 (p(x,y)) 做GMM建模。预测时给定一个新的输入 (x^)从联合分布里取出条件分布 (p(y|x^))。这样做的好处是模型不是全局一条直线而是多个局部线性模型的加权组合。每个高斯分量相当于一个局部工作点在这个工作点附近(x) 和 (y) 近似线性离工作点远了这个分量的权重自动变小预测由其他更相关的分量接管。4.2 条件分布怎么算把第 (k) 个分量的均值与协方差分块[ \mu_k\begin{bmatrix}\mu_k^x\\mu_k^y\end{bmatrix},\quad \Sigma_k\begin{bmatrix}\Sigma_k^{xx} \Sigma_k^{xy}\\Sigma_k^{yx} \Sigma_k^{yy}\end{bmatrix} ]那么在分量 (k) 下给定 (x^*)(y) 的条件分布仍是高斯分布[ \mu_k^{y|x}\mu_k^y\Sigma_k^{yx}(\Sigma_k^{xx})^{-1}(x^*-\mu_k^x) ][ \Sigma_k^{y|x}\Sigma_k^{yy}-\Sigma_k^{yx}(\Sigma_k^{xx})^{-1}\Sigma_k^{xy} ]这和多元高斯条件分布公式一致。每个分量给出一个线性的条件均值斜率就是 (\Sigma_k^{yx}(\Sigma_k^{xx})^{-1})。对于一个新的输入 (x^*)还要计算每个分量的混合权重但这里不能用固定的 (\pi_k)因为不同输入靠近哪个分量不一样。要把 (\pi_k) 乘以输入侧的高斯密度再归一化[ \beta_k(x^)\frac{\pi_k,\mathcal{N}(x^|\mu_k^x,\Sigma_k^{xx})}{\sum_j\pi_j,\mathcal{N}(x^*|\mu_j^x,\Sigma_j^{xx})} ]最终条件均值是局部均值的加权和[ \hat{y}\sum_{k1}^{K}\beta_k(x^*)\mu_k^{y|x} ]条件方差则要同时考虑分量内部方差和分量之间的均值偏移[ \mathrm{Var}(y|x^)\sum_{k1}^{K}\beta_k(x^)\left(\Sigma_k^{y|x}(\mu_k^{y|x}-\hat{y})^2\right) ]这个公式在做置信区间时很有用。很多教程只给条件均值不给方差实际工程里方差不光能画误差带还能用来做自适应控制里的置信度判断。4.3 完整MATLAB示例分段线性数据拟合先造一段分段线性数据来验证GMRrng(1); x (0:0.05:10); y zeros(length(x), 1); idx1 x 4; idx2 x 4 x 7; idx3 x 7; y(idx1) 0.5 * x(idx1) 0.5 0.1 * randn(sum(idx1), 1); y(idx2) 1.5 * x(idx2) - 3.5 0.15 * randn(sum(idx2), 1); y(idx3) 0.8 * x(idx3) 1.4 0.1 * randn(sum(idx3), 1); X [x, y]; K 3; gm fitgmdist(X, K, RegularizationValue, 1e-5, ... Options, statset(MaxIter, 500));预测函数可以单独保存成gmr_predict.mfunction [ymu, yvar, beta] gmr_predict(gm, x_new) % 输入gm为gmdistribution对象x_new为n行q列 % 假设训练数据最后一维是y前面是x q size(x_new, 2); D size(gm.mu, 2); K gm.NumComponents; n size(x_new, 1); ymu zeros(n, 1); yvar zeros(n, 1); beta zeros(n, K); for i 1:n mu_yx zeros(K, 1); var_yx zeros(K, 1); px zeros(K, 1); for k 1:K mux gm.mu(k, 1:q); muy gm.mu(k, q1:D); S gm.Sigma(:, :, k); Sxx S(1:q, 1:q); Sxy S(1:q, q1:D); Syx S(q1:D, 1:q); Syy S(q1:D, q1:D); mu_yx(k) muy Syx / Sxx * (x_new(i, :) - mux); var_yx(k) Syy - Syx / Sxx * Sxy; px(k) gm.ComponentProportion(k) * mvnpdf(x_new(i, :), mux, Sxx); end beta(i, :) px / sum(px); ymu(i) beta(i, :) * mu_yx; % 混合条件方差 yvar(i) 0; for k 1:K yvar(i) yvar(i) beta(i, k) * (var_yx(k) (mu_yx(k) - ymu(i))^2); end end end调用并可视化x_new (0:0.1:10); [ymu, yvar, beta] gmr_predict(gm, x_new); figure; plot(x, y, .); hold on; plot(x_new, ymu, r-, LineWidth, 2); fill([x_new; flipud(x_new)], ... [ymu 1.96 * sqrt(yvar); flipud(ymu - 1.96 * sqrt(yvar))], ... r, FaceAlpha, 0.2, EdgeColor, none); xlabel(x); ylabel(y); legend(训练数据, GMR均值, 95%置信区间, Location, northwest);跑完会发现GMR在三段斜率不同的位置自动切换中间过渡区也不会有明显突变。这就是GMR的优势它把分段线性拟合问题转化成了联合密度估计问题不需要显式知道拐点在哪里也不需要手动设计切换逻辑。5. 常见问题与排错技巧实录5.1 协方差矩阵奇异或接近奇异这是用GMM最常碰到的错误。典型提示是协方差矩阵必须为正定矩阵或者矩阵接近奇异。原因通常是某个分量的样本太少、数据维度太高、存在强相关特征或者K设置得过大。解决办法有几条调大RegularizationValue比如从1e-6调到1e-4相当于在每个分量协方差对角线上加一个小量。换成CovarianceType, diagonal省去特征间协方差稳定性会好很多。对数据做标准化X zscore(X);让每个特征尺度一致。检查K是否明显超过数据里的自然簇数。我在手写EM时也遇到过协方差变奇异的问题所以在M步里给Sigma加一个1e-6 * eye(D)的保底项这比反复调试崩溃体验好得多。5.2 分量数K怎么选GMM的分量数不是越多越好。K太大会过拟合K太小拟合不了多峰结构。我一般用BIC或AIC选KmaxK 10; bic zeros(maxK, 1); aic zeros(maxK, 1); for k 1:maxK gm_k fitgmdist(X, k, RegularizationValue, 1e-5); n size(X, 1); bic(k) -2 * gm_k.NegativeLogLikelihood gm_k.NumParameters * log(n); aic(k) -2 * gm_k.NegativeLogLikelihood 2 * gm_k.NumParameters; end [~, k_bic] min(bic); [~, k_aic] min(aic);画出来看折线图找拐点或者最小值。BIC对模型复杂度惩罚更重实际使用中通常比AIC选出来的K更小更不容易过拟合。5.3 收敛慢、陷入局部最优GMM的EM对初值很敏感。解决办法fitgmdist支持Start参数默认是用K-means初始化比较可靠如果想更稳妥设置Options里的TolFun调小。手写EM时多跑几个随机种子选择对数似然最高的结果。可以在EM迭代过程中增加MaxIter但很多时候局部最优不是迭代不够而是初始化不好K-means基本够用。5.4 数据行列方向放反这是MATLAB里出现频率极高的低端错误。GMM要求每一行是一个样本每一列是一个特征。如果习惯把特征排成行fitgmdist依然会运行但结果完全错乱样本数变成维度协方差维度变成几万乘几万直接内存爆掉。做任何建模之前先size(X)确认第一个维度是样本数第二个维度是特征数。5.5 GMR外推不可靠高斯混合回归的每个分量本质上是局部线性模型拟合范围只在训练数据的输入区间内有效。一旦x_new超出训练范围某个分量的权重可能被推向极端条件均值会沿着该分量的斜率无约束外推输出值很容易发疯。所以实际使用中尽量只预测训练数据覆盖范围内的输入如果必须外推至少把yvar的估计结果也输出用置信区间宽度提醒自己结果不可信。6. 实操心得把这套流程固化成自己的工具箱我最早接触GMM是在做系统辨识项目时需要把非线性对象划分成多个工作点每个工作点建一个局部线性模型。用GMM对输入输出的联合数据进行聚类之后后验概率天然就是这个样本在各个工作点的软权重比人工切分工况段平滑得多。配合前面提到的GMR还能做到不需要先分类、直接对局部模型加权输出整个流程非常自然。有几点个人经验值得分享。第一如果数据量很大注意mvnpdf在循环里容易成为性能瓶颈尽量用fitgmdist的向量化实现手写EM只适合学习和小规模验证。第二GMM出来的聚类标签不能直接当作真实类别它只能反映数据分布的簇结构业务含义需要自己解释。第三做GMR前对输入输出都做标准化能在数值稳定性上省很多心但预测结果要记得反标准化回去。高斯混合模型和高斯混合回归的组合在轨迹预测、机器人运动建模、时间序列分段、系统辨识等方向都有应用场景。把这套MATLAB代码整理成函数训练和预测分开换个数据集就能直接复用也算是一个踏实的工程积累。本文还有配套的精品资源点击获取
高斯混合模型与高斯混合回归：EM算法及MATLAB实现详解
简介高斯混合模型GMM与高斯混合回归GMR是机器学习中常用的概率建模工具前者用于无监督聚类和密度估计后者在此基础上实现连续变量的回归预测。本资源提供一套完整的MATLAB实现代码目标使用者是具备基础统计知识、正在学习机器学习或需要快速搭建GMM/GMR模型的研究人员和学生。压缩包内共26个文件其中19个.m脚本覆盖了高斯概率计算、EM算法中的E步与M步、K均值初始化、混合回归训练与预测、LQR控制等核心流程6个.mat数据文件用于直接运行和复现演示另有1个txt说明文档便于理清代码结构。资源包大小仅98KB轻量便携。目前已有1958人学习代码中配有实例和可视化图形通过运行demo可直观理解参数估计迭代过程和回归效果。借助该资源可以系统掌握GMM聚类的数学原理、GMR回归的预测逻辑以及MATLAB工程化实现思路是兼顾理论学习与实际建模的优质参考。 做数据分析的同行应该都有过这种经历数据分布明明有好几个山头硬拿一个高斯分布去拟合估计出来的均值落在两座山之间的谷底方差也大得离谱完全不能代表真实结构。这时候就得搬出高斯混合模型GMM。用K个高斯分布加权叠加去逼近任意复杂的多峰密度聚类、异常检测、密度估计、回归预测量它都能干。而高斯混合回归则是把GMM从描述数据升级成预测输出在软切换的局部线性回归任务里比单纯套一条全局直线要实用得多。这篇文章不绕弯子直接从原理讲到MATLAB实现把GMM、EM迭代、高斯混合回归的预测函数、常见坑位一次说清楚适合正在做模式识别、信号处理、系统辨识和机器学习相关工作的同学参考。1. GMM到底在算什么从单峰到多峰的概率视角1.1 一个高斯分布为什么不够用高斯分布是统计建模里最常用的假设均值和协方差两个东西就能完整描述一个对称的钟形分布。但实际问题里数据往往不是一个整体。比如同一个工艺线上生产出来的零件可能来自两台精度不同的机床传感器采集的振动信号可能包含正常工况和异常工况两种状态。把这些数据混在一起看直方图大概率会出现两个甚至多个峰。如果用单一高斯分布去拟合这种多峰数据最大似然估计出来的均值会把所有峰值揉在一起结果是既不在这也不在那。协方差也会被迫拉大覆盖两个峰之间的空白区域。密度估计严重失真后续的检测、分类都会跟着出错。这时候把数据拆成几个高斯分量的叠加让每个分量负责描述一个局部聚集结构再按权重组合起来就能得到非常接近真实分布的密度函数——这就是GMM存在的意义。1.2 高斯混合模型的数学表达GMM的密度函数形式上很干净[ p(x)\sum_{k1}^{K}\pi_k,\mathcal{N}(x|\mu_k,\Sigma_k) ]其中 (K) 是分量个数(\pi_k) 是第 (k) 个分量的混合系数满足 (\pi_k\ge 0) 且 (\sum_k\pi_k1)(\mu_k) 和 (\Sigma_k) 分别是第 (k) 个高斯分量的均值向量和协方差矩阵。可以把 (\pi_k) 理解为数据里大约有 (\pi_k) 比例的点是由第 (k) 个分量产生的。举个最简单的例子两个分量(\pi_10.6)(\pi_20.4)第一条高斯中心在 ((1,1))第二条中心在 ((5,5))那么最终生成的数据会呈现两个椭圆簇一个密一些一个疏一些。这里的隐藏变量就是每个样本到底属于哪个分量这个变量我们观测不到所以才需要EM算法来迭代求解。1.3 什么时候优先选GMMGMM不是万能的但它有几个特别合适的场景密度估计面对多峰分布需要给任意一个点输出概率密度值。软聚类像K-means那种硬分配每个点只能属于一个簇GMM给出的是属于每个簇的后验概率对重叠簇更友好。异常检测用训练好的GMM计算样本的似然概率概率低于阈值就判为异常。作为回归或降维的前置模块高斯混合回归、高斯混合隐变量模型等都是在GMM的基本框架上加任务。如果只是做大规模硬聚类K-means更快更省事但当你需要概率、需要椭圆形的簇边界、需要把聚类结果进一步接回归或分类时GMM是更合适的底座。2. EM算法GMM参数估计的核心引擎2.1 为什么极大似然直接做不动GMM的参数包括 (\pi_k,\mu_k,\Sigma_k)。如果直接用极大似然估计需要对数似然函数[ \log p(X)\sum_{i1}^{N}\log\left(\sum_{k1}^{K}\pi_k,\mathcal{N}(x_i|\mu_k,\Sigma_k)\right) ]问题就出在 (\log) 里面套了一个求和。求导之后混合系数的更新会和所有未知参数纠缠在一起得不到闭式解。于是EM算法出场先假设参数已知估计每个样本的隐变量再把隐变量当作已知重新估计参数。两步交替不断抬升对数似然。2.2 E步算责任值在第 (t) 次迭代中用当前参数 ({\pi_k^{(t)},\mu_k^{(t)},\Sigma_k^{(t)}}) 计算第 (i) 个样本属于第 (k) 个分量的后验概率也叫责任值[ \gamma_{ik}\frac{\pi_k^{(t)},\mathcal{N}(x_i|\mu_k^{(t)},\Sigma_k^{(t)})}{\sum_{j1}^{K}\pi_j^{(t)},\mathcal{N}(x_i|\mu_j^{(t)},\Sigma_j^{(t)})} ]这个式子的分母是所有分量的加权密度之和归一化后(\gamma_{ik}) 表示样本 (x_i) 有多大概率来自分量 (k)。如果在某个分量附近责任值就趋近1如果落在两个分量中间责任值会按两边的密度占比分配。2.3 M步按责任值加权更新参数拿到责任值后参数更新变得很直观[ \pi_k^{(t1)}\frac{1}{N}\sum_{i1}^{N}\gamma_{ik} ][ \mu_k^{(t1)}\frac{\sum_{i1}^{N}\gamma_{ik}x_i}{\sum_{i1}^{N}\gamma_{ik}} ][ \Sigma_k^{(t1)}\frac{\sum_{i1}^{N}\gamma_{ik}(x_i-\mu_k^{(t1)})(x_i-\mu_k^{(t1)})^T}{\sum_{i1}^{N}\gamma_{ik}} ]本质上就是加权平均。责任值大的样本对当前分量的均值和协方差贡献大责任值小的贡献小。迭代交替进行直到对数似然的增量小于某个阈值或者达到最大迭代次数。我手写EM时初始化一般用K-means的聚类结果把K-means得到的簇心作为 (\mu) 的初值按簇内样本协方差作为 (\Sigma) 初值(\pi) 设为各类样本占比。这样比随机初始化收敛稳定得多也能有效避开一部分局部最优。3. MATLAB实现GMM现成函数与手写EM两条路3.1 最省事的方式fitgmdistMATLAB里做GMM第一选择是统计工具箱自带的fitgmdist。它内部已经把初始化、EM迭代、协方差类型选择、正则化都封装好了。基本用法rng(42); % 生成两个簇的数据 X [randn(300, 2) * 0.6 [1 1]; randn(200, 2) * 1.2 [5 5]]; k 2; gm fitgmdist(X, k, ... CovarianceType, full, ... RegularizationValue, 1e-6);fitgmdist返回的gm是gmdistribution对象直接用pdf(gm, X)算密度用cluster(gm, X)做硬聚类用random(gm, n)生成服从该混合模型的样本。CovarianceType可以选full或diagonal前者能拟合椭圆形的簇后者计算更稳更快RegularizationValue会在协方差矩阵对角线上加一个小常数防止矩阵奇异这个参数在样本少、维度高时特别关键。3.2 手写一个简化版EM只调现成函数容易踩到不知道内部发生了什么的坑。下面这段代码是简化版EM适合学习function [mu, Sigma, pi_, llh] my_gmm_em(X, K, maxIter, tol) % 输入X为N行D列的数据K为分量数 % 输出mu为K行D列Sigma为D×D×K数组pi_为K维向量 N size(X, 1); D size(X, 2); % 用kmeans做初始化 [idx, C] kmeans(X, K, Replicates, 5); mu C; pi_ zeros(K, 1); Sigma zeros(D, D, K); for k 1:K Xk X(idx k, :); pi_(k) size(Xk, 1) / N; Sigma(:, :, k) cov(Xk) 1e-6 * eye(D); end llh zeros(maxIter, 1); for iter 1:maxIter % E步 resp zeros(N, K); for k 1:K resp(:, k) pi_(k) * mvnpdf(X, mu(k, :), Sigma(:, :, k)); end sum_resp sum(resp, 2); llh(iter) sum(log(sum_resp)); resp resp ./ sum_resp; % 归一化成责任值 % 检查收敛 if iter 1 abs(llh(iter) - llh(iter-1)) tol llh llh(1:iter); break; end % M步 Nk sum(resp, 1); % 1×K pi_ Nk / N; mu (resp * X) ./ Nk; % K×D for k 1:K Xc X - mu(k, :); Sigma(:, :, k) (Xc .* resp(:, k)) * Xc / Nk(k) 1e-6 * eye(D); end end end这个版本的M步里Xc .* resp(:, k)用到了MATLAB的隐式扩展Xc是 (N\times D) 的零均值样本矩阵resp(:, k)是 (N\times 1) 的权重列向量乘出来之后每一行都乘上对应样本的责任值。再左乘Xc就是在计算加权外积和最后除以 (N_k)。如果你用旧版MATLAB可能不支持隐式扩展可以改成repmat。3.3 聚类结果可视化硬聚类只用一行idx cluster(gm, X); gscatter(X(:, 1), X(:, 2), idx);如果想画每个分量的置信椭圆需要从gm.mu和gm.Sigma里取均值与协方差再在椭圆方程里算角度和半轴长度。一个常用思路是对协方差矩阵做特征值分解特征值开方后乘以对应方向的单位向量就得到椭圆的长短轴。配合gscatter出来的图会专业很多。4. 高斯混合回归从联合分布到条件预测4.1 核心思路把y也当成高斯分量的一部分高斯混合回归Gaussian Mixture Regression简称GMR和普通回归的建模路径不一样。普通回归直接建模 (p(y|x))而GMR先把输入 (x) 和输出 (y) 拼成一个向量对联合分布 (p(x,y)) 做GMM建模。预测时给定一个新的输入 (x^)从联合分布里取出条件分布 (p(y|x^))。这样做的好处是模型不是全局一条直线而是多个局部线性模型的加权组合。每个高斯分量相当于一个局部工作点在这个工作点附近(x) 和 (y) 近似线性离工作点远了这个分量的权重自动变小预测由其他更相关的分量接管。4.2 条件分布怎么算把第 (k) 个分量的均值与协方差分块[ \mu_k\begin{bmatrix}\mu_k^x\\mu_k^y\end{bmatrix},\quad \Sigma_k\begin{bmatrix}\Sigma_k^{xx} \Sigma_k^{xy}\\Sigma_k^{yx} \Sigma_k^{yy}\end{bmatrix} ]那么在分量 (k) 下给定 (x^*)(y) 的条件分布仍是高斯分布[ \mu_k^{y|x}\mu_k^y\Sigma_k^{yx}(\Sigma_k^{xx})^{-1}(x^*-\mu_k^x) ][ \Sigma_k^{y|x}\Sigma_k^{yy}-\Sigma_k^{yx}(\Sigma_k^{xx})^{-1}\Sigma_k^{xy} ]这和多元高斯条件分布公式一致。每个分量给出一个线性的条件均值斜率就是 (\Sigma_k^{yx}(\Sigma_k^{xx})^{-1})。对于一个新的输入 (x^*)还要计算每个分量的混合权重但这里不能用固定的 (\pi_k)因为不同输入靠近哪个分量不一样。要把 (\pi_k) 乘以输入侧的高斯密度再归一化[ \beta_k(x^)\frac{\pi_k,\mathcal{N}(x^|\mu_k^x,\Sigma_k^{xx})}{\sum_j\pi_j,\mathcal{N}(x^*|\mu_j^x,\Sigma_j^{xx})} ]最终条件均值是局部均值的加权和[ \hat{y}\sum_{k1}^{K}\beta_k(x^*)\mu_k^{y|x} ]条件方差则要同时考虑分量内部方差和分量之间的均值偏移[ \mathrm{Var}(y|x^)\sum_{k1}^{K}\beta_k(x^)\left(\Sigma_k^{y|x}(\mu_k^{y|x}-\hat{y})^2\right) ]这个公式在做置信区间时很有用。很多教程只给条件均值不给方差实际工程里方差不光能画误差带还能用来做自适应控制里的置信度判断。4.3 完整MATLAB示例分段线性数据拟合先造一段分段线性数据来验证GMRrng(1); x (0:0.05:10); y zeros(length(x), 1); idx1 x 4; idx2 x 4 x 7; idx3 x 7; y(idx1) 0.5 * x(idx1) 0.5 0.1 * randn(sum(idx1), 1); y(idx2) 1.5 * x(idx2) - 3.5 0.15 * randn(sum(idx2), 1); y(idx3) 0.8 * x(idx3) 1.4 0.1 * randn(sum(idx3), 1); X [x, y]; K 3; gm fitgmdist(X, K, RegularizationValue, 1e-5, ... Options, statset(MaxIter, 500));预测函数可以单独保存成gmr_predict.mfunction [ymu, yvar, beta] gmr_predict(gm, x_new) % 输入gm为gmdistribution对象x_new为n行q列 % 假设训练数据最后一维是y前面是x q size(x_new, 2); D size(gm.mu, 2); K gm.NumComponents; n size(x_new, 1); ymu zeros(n, 1); yvar zeros(n, 1); beta zeros(n, K); for i 1:n mu_yx zeros(K, 1); var_yx zeros(K, 1); px zeros(K, 1); for k 1:K mux gm.mu(k, 1:q); muy gm.mu(k, q1:D); S gm.Sigma(:, :, k); Sxx S(1:q, 1:q); Sxy S(1:q, q1:D); Syx S(q1:D, 1:q); Syy S(q1:D, q1:D); mu_yx(k) muy Syx / Sxx * (x_new(i, :) - mux); var_yx(k) Syy - Syx / Sxx * Sxy; px(k) gm.ComponentProportion(k) * mvnpdf(x_new(i, :), mux, Sxx); end beta(i, :) px / sum(px); ymu(i) beta(i, :) * mu_yx; % 混合条件方差 yvar(i) 0; for k 1:K yvar(i) yvar(i) beta(i, k) * (var_yx(k) (mu_yx(k) - ymu(i))^2); end end end调用并可视化x_new (0:0.1:10); [ymu, yvar, beta] gmr_predict(gm, x_new); figure; plot(x, y, .); hold on; plot(x_new, ymu, r-, LineWidth, 2); fill([x_new; flipud(x_new)], ... [ymu 1.96 * sqrt(yvar); flipud(ymu - 1.96 * sqrt(yvar))], ... r, FaceAlpha, 0.2, EdgeColor, none); xlabel(x); ylabel(y); legend(训练数据, GMR均值, 95%置信区间, Location, northwest);跑完会发现GMR在三段斜率不同的位置自动切换中间过渡区也不会有明显突变。这就是GMR的优势它把分段线性拟合问题转化成了联合密度估计问题不需要显式知道拐点在哪里也不需要手动设计切换逻辑。5. 常见问题与排错技巧实录5.1 协方差矩阵奇异或接近奇异这是用GMM最常碰到的错误。典型提示是协方差矩阵必须为正定矩阵或者矩阵接近奇异。原因通常是某个分量的样本太少、数据维度太高、存在强相关特征或者K设置得过大。解决办法有几条调大RegularizationValue比如从1e-6调到1e-4相当于在每个分量协方差对角线上加一个小量。换成CovarianceType, diagonal省去特征间协方差稳定性会好很多。对数据做标准化X zscore(X);让每个特征尺度一致。检查K是否明显超过数据里的自然簇数。我在手写EM时也遇到过协方差变奇异的问题所以在M步里给Sigma加一个1e-6 * eye(D)的保底项这比反复调试崩溃体验好得多。5.2 分量数K怎么选GMM的分量数不是越多越好。K太大会过拟合K太小拟合不了多峰结构。我一般用BIC或AIC选KmaxK 10; bic zeros(maxK, 1); aic zeros(maxK, 1); for k 1:maxK gm_k fitgmdist(X, k, RegularizationValue, 1e-5); n size(X, 1); bic(k) -2 * gm_k.NegativeLogLikelihood gm_k.NumParameters * log(n); aic(k) -2 * gm_k.NegativeLogLikelihood 2 * gm_k.NumParameters; end [~, k_bic] min(bic); [~, k_aic] min(aic);画出来看折线图找拐点或者最小值。BIC对模型复杂度惩罚更重实际使用中通常比AIC选出来的K更小更不容易过拟合。5.3 收敛慢、陷入局部最优GMM的EM对初值很敏感。解决办法fitgmdist支持Start参数默认是用K-means初始化比较可靠如果想更稳妥设置Options里的TolFun调小。手写EM时多跑几个随机种子选择对数似然最高的结果。可以在EM迭代过程中增加MaxIter但很多时候局部最优不是迭代不够而是初始化不好K-means基本够用。5.4 数据行列方向放反这是MATLAB里出现频率极高的低端错误。GMM要求每一行是一个样本每一列是一个特征。如果习惯把特征排成行fitgmdist依然会运行但结果完全错乱样本数变成维度协方差维度变成几万乘几万直接内存爆掉。做任何建模之前先size(X)确认第一个维度是样本数第二个维度是特征数。5.5 GMR外推不可靠高斯混合回归的每个分量本质上是局部线性模型拟合范围只在训练数据的输入区间内有效。一旦x_new超出训练范围某个分量的权重可能被推向极端条件均值会沿着该分量的斜率无约束外推输出值很容易发疯。所以实际使用中尽量只预测训练数据覆盖范围内的输入如果必须外推至少把yvar的估计结果也输出用置信区间宽度提醒自己结果不可信。6. 实操心得把这套流程固化成自己的工具箱我最早接触GMM是在做系统辨识项目时需要把非线性对象划分成多个工作点每个工作点建一个局部线性模型。用GMM对输入输出的联合数据进行聚类之后后验概率天然就是这个样本在各个工作点的软权重比人工切分工况段平滑得多。配合前面提到的GMR还能做到不需要先分类、直接对局部模型加权输出整个流程非常自然。有几点个人经验值得分享。第一如果数据量很大注意mvnpdf在循环里容易成为性能瓶颈尽量用fitgmdist的向量化实现手写EM只适合学习和小规模验证。第二GMM出来的聚类标签不能直接当作真实类别它只能反映数据分布的簇结构业务含义需要自己解释。第三做GMR前对输入输出都做标准化能在数值稳定性上省很多心但预测结果要记得反标准化回去。高斯混合模型和高斯混合回归的组合在轨迹预测、机器人运动建模、时间序列分段、系统辨识等方向都有应用场景。把这套MATLAB代码整理成函数训练和预测分开换个数据集就能直接复用也算是一个踏实的工程积累。本文还有配套的精品资源点击获取