Matlab边缘分布拟合与KS检验：拟合优度检验全流程实战

边缘分布拟合和KS检验是数据分析里几乎绕不开的两个环节。尤其在金融风险度量、水文频率分析、可靠性寿命分析这类场景里拿到观测数据的第一件事往往不是急着上模型而是先把每个变量的边缘分布彻底搞清楚。我最近在项目里用Matlab把“分布拟合 KS检验”这条流程完整跑了一遍从fitdist参数估计到kstest拟合优度检验再到多个候选分布之间的对比选优中间踩了不少坑也整理出一些能直接参考的经验。这篇就围绕边缘分布拟合和KS检验这两个核心点把Matlab里的实现路径、统计原理、常见翻车点一次讲透。适合正在做统计建模、仿真数据生成或者准备进入Copula建模的朋友参考。1. 边缘分布拟合是什么为什么它是建模的第一步1.1 边缘分布与联合分布的关系边缘分布这个词在多元统计分析里非常常见。如果你有一组二维或更高维的随机向量比如同时观测了某流域的日降雨量和洪峰流量那么每个分量单独拿出来看它的概率分布就叫边缘分布。之所以叫“边缘”是因为联合分布函数把多个变量放在一起刻画而单独看某一个变量时相当于把联合概率密度对其他变量做积分积分的结果落在联合分布表的“边缘”上于是就有了这个叫法。在很多实际分析场景里我们一开始并不需要完整的联合分布。主流做法是先分别拟合每个变量的边缘分布再通过Copula或者其他方法把变量之间的相依结构刻画出来。这也是Copula建模的标准流程先定边缘再定相关结构。边缘分布拟合得准不准直接影响后面所有环节的可信度。1.2 拟合的本质找“最像”的那个分布族成员分布拟合的本质是用一个形式已知的概率分布正态、对数正态、Weibull、Gamma等去近似观测数据的真实分布。这里的“真实分布”永远是未知的我们能做的只是在这个已知分布族里面找到最接近真实情况的那一个。这个过程包含两个目标一是把参数估准二是对“拟合得好不好”给出一个可量化的判断。可量化判断这一步就是拟合优度检验要做的事。KS检验Kolmogorov-Smirnov检验是其中用得最广的一种它通过比较经验分布函数和理论分布函数之间的距离来判断样本是否来自某个指定分布。参数估计决定“我们拟合的是哪一条曲线”KS检验决定“这条曲线到底像不像数据”。1.3 边缘分布拟合的典型应用场景Copula建模的前置步骤金融、水文、气象等领域常用Copula刻画变量间的相依结构。每个变量的边缘分布先单独拟合好再做概率积分变换得到均匀分布序列才能进入下一步的Copula参数估计。蒙特卡洛模拟做风险分析或者系统仿真时需要从一个已知分布里抽样。没有拟合好的边缘分布模拟就无从谈起。预测区间构造回归和机器学习模型的残差经常需要假设服从某个分布。残差的边缘分布是否成立直接决定预测区间的可靠性。可靠性寿命分析设备寿命数据通常用Weibull、Gamma等分布拟合用来估算失效率、MTBF等关键指标。不管哪种场景套路都是一样的先把数据拟合成分布再检验拟合效果。所以这篇文章的重点就放在“拟合 检验”这两个动作上。2. Matlab分布拟合工具箱怎么用fitdist与mle2.1 fitdist一行代码完成MLE参数估计Matlab里做分布拟合绝大多数情况用fitdist就够了。它会用极大似然估计MLE自动估计参数并且返回一个概率分布对象。这个对象不是简单的一堆数字它把参数、标准误差、AIC、BIC等信息都打包在了一起。% 拟合正态分布 pd fitdist(x, Normal); % 查看参数和拟合评价指标 mu_hat pd.mu; sigma_hat pd.sigma; aic_val pd.AIC; bic_val pd.BIC;我最喜欢fitdist的一点是返回的分布对象可以直接被其他函数接收。比如cdf(pd, x)可以算理论累积概率random(pd, n, 1)可以按拟合结果生成随机样本qqplot(x, pd)可以直接画Q-Q图。这种面向对象的接口设计让整个工作流非常顺。如果不知道怎么选分布还有一个很实用的GUI工具在命令行敲dfittool老版本叫Distribution Fitting Tool新版本通过App标签页里的Distribution Fitter打开。它能交互式地让你点选不同分布族即时查看拟合曲线和统计量。我建议新手先用这个工具做一次探索性分析心里有数了再写脚本批量跑。2.2 常用分布速查表不同分布的适用场景差异很大选错分布是新手最容易犯的错误。我整理了一张速查表方便对照选择分布名称fitdist中的名称参数典型适用场景正态分布Normalmu, sigma对称数据、回归残差、测量误差对数正态Lognormalmu, sigma正数且右偏的数据、金融资产价格WeibullWeibullA尺度, B形状设备寿命、可靠性、风速数据GammaGammaa形状, b尺度降雨量、等待时间、索赔金额指数分布Exponentialmu均值无记忆性事件间隔、排队论极值分布ExtremeValuemu, sigma极值分析、环境监测峰值t位置尺度tLocationScalemu, sigma, nu厚尾数据、金融收益率均匀分布Uniformlower, upper圆形或截断数据、无信息先验对数正态和Gamma都适合正数右偏数据但两者在尾部的表现差异很大Fisher准则下不容易区分常需要借助业务背景判断。后面第4章我会演示怎么用统计数据帮忙做选择但最终拍板时还是要结合数据来源的物理意义。2.3 mle函数遇到自定义分布时的替代方案fitdist覆盖了Matlab内置的三十多种分布绝大多数场景够用了。但如果你需要拟合一个内置列表里不存在的分布比如带位置参数的分布、混合分布或者某个行业专用的分布族就得请出mle函数。mle可以直接接收自定义的概率密度函数句柄也支持指定初值和边界。% 用mle拟合自定义Burr分布Burr XII型 pdf_burr (x, a, c, k) ... (c * k ./ a) .* (x ./ a).^(c - 1) ./ (1 (x ./ a).^c).^(k 1); params mle(x, pdf, pdf_burr, ... start, [1, 2, 1], ... LowerBound, [1e-6, 1e-6, 1e-6]);这里有几个细节值得注意。第一start初值非常重要MLE的优化过程是迭代求解的初值给得离谱很可能直接不收敛。第二LowerBound和UpperBound可以把参数限制在合理范围内比如尺度参数必须是正数。第三mle返回的是一组参数值而不是分布对象如果需要生成随机数可以用makedist手动建对象或者直接用参数代入random函数。如果你和一样经常分位数那个初值的问题我的习惯是先做一次矩估计把矩估计结果作为MLE的初值收敛速度和稳定性都会好很多。2.4 参数估计中的数值坑数据边界问题Lognormal要求数据严格为正Gamma也一样。如果你的数据里有0拟合时会报错或者给出没有意义的结果。我通常先看数据范围必要时做一个小偏移处理。Weibull拟合不收敛形状参数B的搜索结果对初值比较敏感尤其是样本量小的时候。可以多试几组初值或者用fitdist之前先画个Weibull概率图通过图上散点斜率估算初值。和异常值的爱恨情仇MLE对异常值比较敏感尤其是正态分布。如果数据里有明显的离群点拟合出来的均值和标准差都会被拉偏。我的习惯是先做一次异常值筛查但别机械地删数据要结合业务判断这些点该不该保留。第2章解决了“怎么拟合”的问题接下来要回答“拟合得好不好”也就是KS检验登场的时候。3. KS检验原理与Matlab实现3.1 KS检验的统计思想KS检验的核心思想非常直观比较经验分布函数ECDF和假设的理论分布函数CDF看它们之间的“最大差距”有多大。D sup|F_n(x) - F(x)|其中F_n(x)是样本经验分布F(x)是理论累积分布函数sup表示取上确界也就是两条曲线之间最大的那个垂直距离。如果这个最大距离足够小我们就认为样本来自这个理论分布如果太大就拒绝原假设。KS检验不依赖具体的分布形式不对数据做分组所以它比卡方检验更精细而且在小样本下也能用。这是它最吸引人的地方不需要知道总体分布长什么样只要能把理论CDF算出来就能检验。Matlab里实现KS检验的函数有三个常用接口% 检验数据是否服从标准正态分布 [h, p, ksstat] kstest(x); % 检验数据是否服从指定的分布对象比如拟合好的分布 [h, p, ksstat] kstest(x, CDF, pd); % 检验数据是否服从指定CDF函数 cdf_fun (x) wblcdf(x, A_hat, B_hat); [h, p, ksstat] kstest(x, CDF, cdf_fun);返回结果中h1表示在显著性水平下拒绝原假设h0表示不能拒绝p是p值ksstat就是刚才公式里的D统计量。3.2 kstest与kstest2单样本和双样本的区别刚才说的都是单样本KS检验用来判断“数据是否来自某个指定分布”。但实际工作中还有个高频需求判断两组独立样本是否来自同一总体。这时候要用kstest2% 两组数据是否同分布 [h, p, ksstat] kstest2(x1, x2);kstest2的原理和单样本版本类似只是把理论CDF换成了第二组样本的经验CDF比较两条经验分布曲线之间的最大距离。它不限制两组样本量必须相同也不要求数据服从什么特定分布是两样本分布比较的利器。举两个实际例子。我在一个对比实验里用kstest2判断两种工艺下产品尺寸分布是否一致省去了大量参数对比的麻烦。另一个场景是数据增强前后的特征分布漂移检测kstest2的p值变化能直观反映数据分布有没有被改动。注意kstest2的默认原假设是两个样本来自同一连续分布。如果你的数据是离散的或者有很多重复值检验结果会偏向保守具体处理办法在第5章详细说。3.3 核心坑点估计参数带来的检验保守性这一节是整个项目里我印象最深的部分也是绝大多数教程都不会讲透的地方。如果你用fitdist拟合出一个分布对象然后直接把它传给kstest做检验从表面看流程很顺畅但从统计角度看这里隐藏了一个微妙的问题KS检验的理论临界值是在“分布参数已知”的前提下推导出来的。可你手里的参数是用同一组数据估出来的。参数估计相当于已经“偷看”了一次数据拟合出的分布会天然地向观测数据靠拢于是经验分布和理论分布之间的最大距离D会被系统性地压低。这导致的直接后果是p值偏高检验更容易“不拒绝原假设”也就是检验功效降低。参数越多这个偏差越明显。正态分布有2个参数Weibull也是2个但如果你用tLocationScale这种3参数的分布保守性会更严重。一句话总结直接用拟合后的分布对象做KS检验结果偏乐观容易把“其实不太像”的分布误判为“可以接受”。3.4 怎么修正Lilliefors检验与蒙特卡洛临界值针对正态分布Matlab提供了修正方案lillietest。这是Lilliefors检验专门处理均值和方差从样本中估计的情况。它的临界值是通过模拟得到的比直接用kstest更可靠。% 正态性检验参数未知从样本估计 [h_lf, p_lf] lillietest(x);但Lilliefors只适用于正态分布。如果你拟合的是Weibull、Gamma或者其他分布没有现成函数可用就需要自己算经验临界值。做法是蒙特卡洛模拟从拟合好的分布里反复抽样对每个模拟样本重新估计参数、重新计算KS统计量最后取这些统计量的高分位数作为经验临界值。% 蒙特卡洛计算经验临界值 rng(2024, twister); n_sim 5000; n_data numel(x); ks_sim zeros(n_sim, 1); parfor i 1:n_sim % 从拟合好的分布中生成同长度样本 x_sim random(pd_best, n_data, 1); % 重新估计参数 pd_sim fitdist(x_sim, dist_best); % 计算KS统计量 [~, ~, ks_stat] kstest(x_sim, CDF, pd_sim); ks_sim(i) ks_stat; end alpha 0.05; crit_value quantile(ks_sim, 1 - alpha); fprintf(经验临界值(alpha0.05): %.4f\n, crit_value);然后拿观测数据的KS统计量和这个经验临界值比较如果观测的ksstat大于临界值就拒绝原假设否则不能拒绝。这套做法虽然代码写起来多一些但对任何分布族都成立是我现在处理拟合优度检验的标准方案。4. 完整实操从观测数据到拟合选优的全流程4.1 数据预处理与描述统计理论讲完进入实战。假设我手里有一份某流域的日降水量观测数据共200个样本准备做边缘分布拟合。先做数据清洗和基础探索% 数据加载与基础清洗 data readtable(观测数据.xlsx); x data.Var1; x x(~isnan(x)); % 剔除缺失值 % 基础描述统计 fprintf(样本量 n %d\n, length(x)); fprintf(均值 %.4f, 标准差 %.4f\n, mean(x), std(x)); fprintf(偏度 %.4f, 峰度 %.4f\n, skewness(x), kurtosis(x));偏度为正且数值不小说明数据明显右偏正态分布大概率不合适。峰度大于3说明比正态分布更尖峰厚尾。这个直觉很重要它能帮我在候选分布清单上做第一步筛选。此时先画个直方图叠加核密度估计直观感受一下分布形态figure(Color, w); histogram(x, Normalization, pdf, FaceAlpha, 0.3); hold on; ksdensity(x); legend(直方图, 核密度估计); grid on;4.2 多候选分布的批量拟合基于偏度信息我把候选分布定为Lognormal、Weibull、Gamma、Exponential外加一个正态分布做对照组。写个循环批量拟合并计算KS检验dist_list {Normal, Lognormal, Weibull, Gamma, Exponential}; T table(); for k 1:numel(dist_list) try pd fitdist(x, dist_list{k}); T.Distribution(k) dist_list{k}; T.AIC(k) pd.AIC; T.BIC(k) pd.BIC; [h, p] kstest(x, CDF, pd); T.KS_p(k) p; T.H(k) h; catch T.Distribution(k) dist_list{k}; T.AIC(k) NaN; T.BIC(k) NaN; T.KS_p(k) NaN; T.H(k) NaN; end end % 按AIC升序排列 T sortrows(T, AIC); disp(T);这里用try...catch包住每个拟合过程避免某个分布因为边界问题直接中断整个循环。实际跑出来的结果大概是这样的分布AICBICKS p值是否拒绝H0Gamma465.9472.60.372否Lognormal468.2474.90.214否Weibull471.5478.20.106否Normal512.7519.40.018是Exponential520.3523.60.001是先看KS p值这一列正态和指数都被拒绝了和预期一致。Gamma、Lognormal、Weibull三个都没被拒绝单靠KS检验暂时分不出高下。这时候就要靠AIC和BIC来帮忙。4.3 图形诊断直方图、Q-Q图、经验CDF统计量只能给出数字结论最终判断还得多看几张图。Q-Q图是最直观的诊断工具如果数据服从目标分布散点应该基本落在直线上。% Q-Q图诊断以Gamma为例 figure(Color, w); qqplot(x, pd_gamma); grid on;经验CDF叠加理论CDF也是我必看的图figure(Color, w); ecdf(x); hold on; x_grid linspace(min(x), max(x), 200); plot(x_grid, cdf(pd_gamma, x_grid), r-, LineWidth, 1.5); legend(经验CDF, Gamma理论CDF, Location, best); grid on;两条曲线贴合得越好说明拟合越可靠。我特别关注尾部区域因为很多业务决策恰恰依赖尾部分布。如果理论CDF在尾部明显偏离经验点哪怕整体KS检验通过这个分布用于尾部风险估计也是不靠谱的。还有一个容易被忽略的操作把Gamma、Lognormal、Weibull三条理论CDF画在一张图里和ECDF对比尾部差异一目了然。这是我选型时最依赖的视觉证据。4.4 综合选优KS p值、AIC/BIC和业务需求选型不能只盯一个指标。KS p值只能说明“不能拒绝”AIC和BIC则可以在多个候选模型之间做相对比较。AIC和BIC都兼顾了拟合优度和参数数量其中BIC对参数数量的惩罚更重。从上面的表格看Gamma的AIC和BIC都是最小的说明它在拟合优度和模型简洁度之间取得了最佳平衡。Lognormal和Weibull的AIC差距不大但BIC差距也小。单看统计量我会选Gamma但还得结合业务背景确认在降水量分析领域Gamma分布是经典选择物理意义也说得通。统计证据加上业务逻辑都指向Gamma这个决定做起来就踏实了。选定Gamma之后还有个收尾动作对最终选定的模型做一次蒙特卡洛修正的KS检验避免第3.4节说的保守性问题。把观测数据的ksstat和模拟出来的经验临界值比较确认结论依然成立这才算真正完成。4.5 边缘分布移交下游概率积分变换边缘分布拟合不只是画几张漂亮的图最终目的是为下游分析提供输入。在Copula建模或者其他多元分析中拟合好的边缘分布需要做一次概率积分变换PIT把原始数据变成均匀分布% 概率积分变换将边缘分布转为均匀分布 u cdf(pd_gamma, x); % 验证变换后的序列是否接近均匀分布 [h_unif, p_unif] kstest(u, CDF, makedist(Uniform, 0, 1));变换后的u变量理论上服从(0,1)均匀分布。这一步的同时也顺带验证了边缘分布拟合的效果如果u的分布明显偏离均匀说明边缘分布拟合有问题需要回去重新检查。这也是我在项目中特别依赖的一个闭环验证方式。5. 常见问题与排查技巧实录5.1 参数估计不收敛或者警告频出现象是fitdist报错“未找到可行的起始值”或者mle迭代超过了最大次数。我排查时先做三件事第一检查数据里是否有NaN、Inf、0或者负数。Lognormal和Gamma对非正数零容忍数据里混入一个0就可能让似然函数崩掉。第二检查初值是否合理。fitdist内部有默认初值策略但遇到极端偏态数据时也可能不给力这时候用mle手动指定初值。第三数据范围是否过大导致数值溢出。降水量数据经常跨越多个数量级对数变换一下或者调整参数尺度可以缓解。% 手动指定初值的稳健做法 pd fitdist(x, Gamma); % 如果失败改用矩估计结果作为初值 mu_x mean(x); var_x var(x); a0 mu_x^2 / var_x; b0 var_x / mu_x; params mle(x, pdf, (x, a, b) gampdf(x, a, b), ... start, [a0, b0]);5.2 所有候选分布都被KS拒绝了这种情况我也遇到过几次。先别急着给数据强行套分布而是想想是不是数据本身就不是“单一分布”能描述的。常见原因有三个第一个是混合分布。数据可能来自两个或多个总体的混合比如设备故障数据里既有早期失效又有磨损失效单峰分布自然拟合不了。此时可以考虑拟合混合模型比如两个Gamma的混合。第二个是数据存在强截断或删失。如果降水量观测有“低于检测限记为0”的情况分布就变成零膨胀分布了。第三个是异常值干扰。回看数据生成的原始过程往往会发现某段时间的观测有明显干扰因素。排查这一步我通常从图形入手画出数据的分组直方图看是否呈现自然的双峰或者多峰结构。如果确认是混合分布那就要换分析框架而不是硬套单分布。5.3 样本量太小怎么办样本量少于30的时候KS检验的功效很低几乎不太可能拒绝原假设但“不拒绝”也说明不了什么。这种情况下我会倾向于用经验分布描述数据不做参数化分布假设如果必须拟合参数分布优先选择参数尽量少的分布族比如指数分布或者2参数的Weibull同时用Bootstrap方法给出参数估计的置信区间别只看点估计图形上多用核密度估计和带状Q-Q图减少小样本的视觉误导。记住一个原则小样本下KS检验的“不拒绝”不等于“证明来自这个分布”它只能说明“数据量和这个分布不算矛盾”。写在分析报告里时措辞要严谨。5.4 数据里有大量重复值带来的结效应KS检验假设连续分布如果数据里出现大量重复值比如记录精度是整数导致很多相同的观测值经验分布函数会出现台阶这会降低D统计量的灵敏度让检验偏向保守。处理办法有两个一是对数据做轻微抖动加上一个远小于观测精度的高斯噪声二是改用针对离散数据的检验方法。抖动操作我一般这样写% 轻微抖动打破结效应 jitter 0.001 * std(x) * randn(size(x)); x_jittered x jitter;抖动之后要记得在方法论部分如实说明保证分析可复现。5.5 问题速查表现象最可能原因处理建议fitdist报错找不到初值数据含非正值或异常偏态清洗数据或手动指定初值KS检验全部通过但图形明显不符样本量太小或参数估计保守用蒙特卡洛临界值重新检验所有分布都被拒绝混合分布、零膨胀或截断换混合模型或重新审视数据来源p值总是很高参数估计导致的保守性改用lillietest或模拟临界值拟合结果不稳定初值敏感或样本量小多组初值对比Bootstrap评估稳定性数据有大量重复值观测精度低或数据离散轻微抖动或换用离散检验6. 实测后的几点心得项目收尾之后复盘有三条经验让我改变了之后所有分析的习惯。第一条永远不要只盯着p值。KS检验的p值只是一个侧面证据尤其是参数从样本里估出来的时候p值天然偏乐观。我现在每做一个边缘分布拟合必定配套看Q-Q图和ECDF叠加图图形不会骗人。第二条拟合得好不等于有业务意义。统计上最优的分布不一定能解释数据的物理背景。选型时我会把统计指标和领域常识放在一起权衡两者冲突时优先怀疑统计指标背后是否有没被发现的数据问题。第三条复现性是整个流程的最后一道保险。所有数据处理步骤、初值选择、抖动参数、随机数种子都要记录在案。我在项目里就因为改过一次抖动幅度导致前后两批结果不可比白白浪费了两天时间才定位到问题。现在我的脚本开头固定写一行rng(2024, twister)所有涉及随机性的环节从根源上锁死。边缘分布拟合和KS检验这套组合拳看起来基础但真正用得扎实并不容易。把参数估计的数值细节、检验的保守性修正、选优的多元证据链都打通之后它带给下游模型的增益是非常可观的。希望这篇能帮你少走一些我走过的弯路。