简介面向需要系统掌握MATLAB高级数理统计功能的科研与工程人员这套压缩包聚焦实战应用内容涉及多变量分析、假设检验、非参数检验、回归与拟合、时间序列分析、随机过程、生存分析、贝叶斯统计以及聚类判别等核心主题可帮助读者突破常规统计建模瓶颈。包内共2个文件含1段mp4视频演示与1个m脚本代码文件视频便于跟随讲解逐步理解统计方法的实现流程脚本可按自己的数据结构直接修改运行从而加深对代码细节的掌握压缩包整体约15.23MB轻量便于下载。目前已有115人学习下载适合希望在数据分析、预测建模与科学计算中提升统计应用能力的读者。文件结构简洁m脚本中的注释与视频讲解相互配合便于理解和复用通过这些内容可同时获得方法讲解、可执行代码与统计思维训练并能快速上手PCA、ARIMA、MCMC等实用技术。1. 从「11 matlab数理统计高级篇.zip」说起一个压缩包背后要补的三块硬功夫拿到一个写着「11 matlab数理统计高级篇.zip」的资源包多数人的第一反应是解压、翻目录、找看起来最像「高级」的那个 .m 文件然后直接 F5 运行。跑不通的占大多数报错信息往往还都很朴素变量未定义、矩阵维度不一致、函数或变量找不到。真正卡住人的从来不是代码本身而是压缩包里默认你已经会的那三块东西——概率分布对象怎么构造和复用、统计检验在不同数据结构下该选哪一个、仿真结果凭什么算可复现。看到 zip 先确认一件事它通常只交付脚本不交付随机种子的约定、不交付参数估计的置信区间口径、也不交付多重比较的校正方法。这个包要能用起来缺的正是这三块硬功夫而不是更多函数名。2. MATLAB 数理统计的分布对象体系与可复现随机数2.1 用 makedist 与 fitdist 把分布变成可传递的对象很多从入门教程过渡到 matlab 数理统计高级篇的人会一直停在normpdf、normcdf、normrnd这类「函数名 分布名前缀」的写法上。这套写法的致命问题是分布一旦换掉整套代码要逐行改前缀分布参数一旦从已知变成待估函数签名也对不上。较新的 MATLAB 版本提供的是面向对象的分布对象体系分布本身变成一个可以赋值、可以传参、可以当结构体字段存下来的变量。% 构造一个参数已知的分布对象后续所有函数复用同一对象 pd makedist(Normal, mu, 0, sigma, 1); % 从样本反推参数fitdist 返回一个已经拟合好的分布对象 x wblrnd(2, 1.5, 2000, 1); % 生成 Weibull 样本A2 B1.5 pdFit fitdist(x, Weibull); % 极大似然估计 disp(pdFit.ParameterNames); % 输出 {A,B}A 尺度B 形状 disp(pdFit.ParameterValues); % 输出估计值 ci paramci(pdFit, Alpha, 0.05); % 参数的 95% 置信区间 disp(ci); % 两列分别对应 A、B 的下上界这段代码做了三件事。makedist负责在参数已知时构造理论分布返回的对象可以直接喂给random、cdf、icdf、pdf也可以作为假设检验里CDF参数的对照。fitdist负责参数未知时的点估计默认走极大似然Distribution位置参数支持Normal、Weibull、Lognormal、Gamma、Kernel等一批名字。paramci是拟合后对象的方法Alpha默认 0.05Parameter还可指定只算某一个参数的区间。提示fitdist在样本被截断时一定要传Censoring否则极大似然会把删失点当成完整观测参数会系统性偏小。Kernel是里面唯一的非参数选项需要额外指定Kernel如epanechnikov、normal和Supportunbounded、positive、unit interval。它的ParameterValues是带宽和核中心不能拿去算 AIC这点在写脚本时最容易踩。2.2 rng 的状态控制与并行环境下的随机源选择数理统计里最容易被忽略、事后最难补救的一环是随机数能不能复现。同一份脚本在别人的机器上跑出不一样的 p 值几乎全部来自随机流没有固定。MATLAB 的rng提供三种典型用法对应三种不同目的。s rng(20240101, twister); % 固定种子并返回当前生成器状态结构体 a randn(1, 5); rng(s); % 把状态还原后续序列完全重放 b randn(1, 5); disp(isequal(a, b)); % 输出 1说明复现成功 rng(shuffle); % 需要每次不同时使用按时间做种子rng第一个参数是种子第二个参数是生成器类型常用取值有twisterMersenne Twister单线程默认、combRecursive组合递归周期极长适合并行、philox和threefry计数器型显式支持并行子流。保存返回的结构体再还原比只记一个种子更可靠因为状态结构体同时记录了生成器和位置。注意在parfor内部直接调用randnMATLAB 会给每个 worker 分配独立的子流以保证统计独立但这意味着结果会随 worker 数量变化。要严格复现先把 worker 数固定再在parfor之前用rng(seed,combRecursive)统一设置主种子。2.3 拟合优度诊断QQ 图、K-S 检验与 AIC 横向比选拿到一组样本后先问「它服从什么分布」这是 matlab 数理统计里出现频率最高的操作。凭直方图目测形状最容易出错尤其是尾部。工程上常用的判断顺序是先画 QQ 图看整体偏离再用 K-S 或 Anderson-Darling 做定量检验最后用 AIC 在若干候选分布之间比选。% 单样本 K-S检验样本是否来自指定理论分布 pd makedist(Normal, mu, mean(x), sigma, std(x)); [h, p, ksstat] kstest(x, CDF, pd); fprintf(K-S h%d p%.4g stat%.4f\n, h, p, ksstat); % 多候选分布的 AIC 比选 names {Normal, Weibull, Lognormal, Gamma}; aic zeros(numel(names), 1); for k 1:numel(names) pdK fitdist(x, names{k}); nll negloglik(pdK); % 负对数似然 aic(k) 2 * nll 2 * numel(pdK.ParameterValues); % AIC 2nll 2k end [~, best] min(aic); fprintf(最优分布%s\n, names{best});kstest的三个返回值分别是拒绝原假设的标志、p 值和检验统计量CDF传分布对象即可对照理论分布h1表示在默认 0.05 水平下拒绝。AIC 部分的关键是negloglik与参数个数的对应关系Normal和Lognormal各两个参数Weibull和Gamma也是两个所以这里自由度一致比的是拟合优度本身。候选分布典型适用场景需要留意的边界Normal测量误差、均值型指标对偏态和厚尾极敏感Weibull寿命、强度、失效时间形状参数小于 1 时尾部很重Lognormal收入、粒径、细胞数取对数后要再做一次正态检验Gamma等待时间、降水量、保险赔付尺度参数解释依赖业务背景提示K-S 检验对参数是用样本估计出来的情况偏保守p 值会偏大。样本量不大时优先用adtest或lillietest也可以直接做参数 bootstrap 的 K-S。3. 假设检验与方差分析从 ttest 到 anovan 的参数怎么设3.1 检验选型一张表定下 ttest2、ranksum 还是 signrank检验函数选错后面结论全部作废而且脚本不会报任何错。选型只看三件事比较的是均值还是分布位置、样本是否配对、方差是否齐性。把这些条件列成表比记函数名有用得多。场景函数调用必须设置的参数单样本均值对常数ttest(x, mu)Tail取both/right/left双样本、方差齐ttest2(x, y)Vartype,equal双样本、方差不齐ttest2(x, y)Vartype,unequalWelch 校正配对样本ttest(x, y)直接传两个等长向量即为配对不满足正态ranksum/signrankmethod取exact或approximate多组独立、非参数kruskalwallisdisplay取on方差不齐时硬用equal在两组样本量差异大的情况下第一类错误率会明显膨胀。稳健做法是先跑一次方差齐性检验再分支代码结构如下。load(carbig.mat, MPG, Origin); us MPG(strcmp(Origin, USA)); jp MPG(strcmp(Origin, Japan)); [~, vp] vartest2(us, jp); % 方差齐性检验只取 p 值 if vp 0.05 vtype unequal; % 方差不齐走 Welch else vtype equal; end [h, p, ci, stats] ttest2(us, jp, Vartype, vtype, Alpha, 0.01); fprintf(h%d t%.3f df%.1f p%.4g CI[%.3f, %.3f]\n, ... h, stats.tstat, stats.df, p, ci(1), ci(2));vartest2的第二个输出才是 p 值第一个是拒绝标志写成[h, vp]再取vp最不容易出错。ttest2返回的stats里有tstat和dfWelch 情形下df是小数这是正常现象。Alpha默认 0.05把它显式写出来是多组比较前控制族错误率的第一步。注意ttest(x, y)与ttest2(x, y)在语法上完全兼容传参错了也不会报错只会给出错误的结论。配对设计必须用ttest独立设计必须用ttest2。3.2 anovan 多因素方差分析与 multcompare 事后比较多因素设计的方差分析用anovan它比anova1、anova2笨一点因为要先把分组变量打包成 cell但换来的好处是支持任意因子数、不平衡设计和指定交互项。% y 为响应向量g1/g2/g3 为分组向量长度必须一致 [p, tbl, stats] anovan(y, {g1, g2, g3}, ... model, interaction, ... % 主效应 所有两两交互 varnames, {温度, 压力, 批次}, ... sstype, 3, ... % 不平衡数据用 III 型平方和 display, on); % 事后两两比较Bonferroni 控制族错误率 [c, m, h, gnames] multcompare(stats, ... CType, bonferroni, ... Alpha, 0.05, ... Display, on);model是这里最关键的参数。linear只做主效应interaction加所有两两交互full加全部交互项也可以直接传一个矩阵来自定义哪些因子参与交互比如[1 0 0; 0 1 0; 0 0 1; 1 1 0]表示三个主效应加上「因子 1 与因子 2」的交互。sstype默认是 3样本量在各组间不均衡时I 型和 II 型平方和会随因子进入顺序变化III 型不会所以不平衡设计务必显式写 3。multcompare的CType决定事后比较的校正方式bonferroni最保守tukey-kramer均衡设计下更紧scheffe适合任意线性组合的比较lsd不做校正只适合事先已经锁定了少数几对比较的情况。返回的c矩阵前两列是要比较的两组编号末列是 p 值m是各组的估计边际均值。提示anovan要求分组变量要么是数值向量要么是 cell 字符串数组。如果原始数据是 table先用categorical转换再取double否则因子水平顺序会按字典序而不是你期望的业务顺序排。3.3 p 值多重校正手写 Benjamini-Hochberg 控制 FDR一次跑几十上百个检验不做校正就报「显著」是统计脚本里最严重的错误。Bonferroni 控制的是族错误率比较次数一多就几乎什么都检不出来工程上更常用的是 Benjamini-Hochberg 控制错误发现率。MATLAB 的统计与机器学习工具箱没有直接提供 BH 函数手写反而更透明。p [0.001, 0.008, 0.039, 0.041, 0.12, 0.33, 0.65]; m numel(p); [sorted_p, idx] sort(p(:)); % 升序排列并记录原位置 bh sorted_p .* m ./ (1:m); % 第 i 个乘 m/i bh flipud(cummin(flipud(bh))); % 从大到小取累计最小保证单调 q zeros(m, 1); q(idx) min(bh, 1); % 按原顺序还原并截断到 1 disp(q);三个步骤对应 BH 的定义。先把 p 值升序排列并记下原始下标第 i 小的 p 值乘以 m/i 得到校正后的 q 值然后从最大的一头往回取累计最小值这一步是为了让 q 值单调不减否则会出现「更小的原始 p 值算出更大的 q 值」这种反直觉结果。最后按原索引还原并截断到 1。判定显著时把 q 与目标 FDR 水平常取 0.05 或 0.10比而不是跟原始 p 值比。如果是基因表达这类场景Bioinformatics Toolbox 里的mafdr直接支持BHFDR选项还额外实现了 Storey 的 q 值方法估计 π0 后检验效力更高。两者结论方向一致但 q 值数值会有差异写报告时要注明用的是哪一种。4. 回归与降维fitlm、fitglm、pca 的工程化用法4.1 fitlm 稳健回归与残差诊断量线性模型在 matlab 数理统计高级篇里出现的频次仅次于假设检验但多数脚本只打印mdl.Coefficients就收工不看残差也不看杠杆点。fitlm的正确用法是把 table 直接传进去让变量名自动进入公式然后用诊断图确认假设是否成立。tbl table(x1, x2, x3, y, VariableNames, {x1,x2,x3,y}); mdl fitlm(tbl, y ~ x1 x2 x3 x1:x2, ... RobustOpts, on, ... % bisquare 稳健拟合抑制离群点 Intercept, true); disp(mdl.Coefficients); % Estimate / SE / tStat / pValue disp(mdl.Rsquared.Adjusted); % 调整 R²比 R² 更可信 anova(mdl); % 各项平方和与 F 检验 plotResiduals(mdl, fitted); % 残差对拟合值看是否喇叭形 plotDiagnostics(mdl, leverage); % 杠杆值找强影响点RobustOpts取on时用 bisquare 权重迭代重加权等价于对离群点降权如果想自定义权重函数传结构体struct(wfun,fair,tune,1.5)即可tune是调优常数。Rsquared是结构体含Ordinary、Adjusted、Ordinary三个字段多元回归里看Adjusted。残差图出现明显的漏斗形说明异方差此时应改用稳健标准误或对响应做变换而不是直接读 p 值。逐步回归用stepwiselmPEnter和PRemove默认分别是 0.05 和 0.10后一个比前一个大是为了防止变量在边界上反复进出。这个函数对样本量小、候选变量多的数据容易过拟合用它筛出来的模型一定要在留出集上重新验证一次。4.2 fitglm 做二分类与链接函数选择响应是 0/1 或者计数时不能再用fitlm改用fitglm。它比旧的glmfit友好在于直接接受 table 和公式并且自带predict和devianceTest。% X 是 n×p 预测矩阵y 是 n×1 的 0/1 向量 mdl fitglm(X, y, Distribution, binomial, ... Link, logit, ... % 也可取 probit、comploglog LikelihoodPenalty, jeffreys-prior); % 完全分离时稳定估计 prob predict(mdl, Xnew); % 输出属于正类的概率 label prob 0.5; % 模型整体显著性与常数模型做偏差检验 devianceTest(mdl);Distribution取binomial时Link可选logit、probit、comploglog取poisson时通常配log。LikelihoodPenalty是较新版本引入的正则化选项jeffreys-prior在样本量小、类别完全可分的情况下能显著稳定系数否则标准误会膨胀到不可用。predict返回的概率不要直接用 0.5 切类别不平衡时应该按业务代价调整阈值或者用 ROC 曲线找约登指数最大的点。如果坚持用glmfit要自己给设计矩阵补一列全 1返回的b是系数、dev是偏差、stats.p是每个系数的 p 值接口更底层但输出更杂。新脚本建议直接用fitglm。4.3 pca 降维与因子分析在特征压缩中的落地维度高了以后回归系数不稳定先降维再建模是常见流程。pca的输出字段多写脚本时按名字取不要依赖位置顺序。[coeff, score, latent, tsquared, explained, mu] pca(X, ... NumComponents, 10, ... % 只保留前 10 个主成分 Algorithm, svd, ... % svd 精度最高eig 更快 Centered, true, ... % 默认按列中心化 VariableWeights, variance); % 变量量纲差异大时按方差加权 cumExp cumsum(explained); k find(cumExp 85, 1); % 取到累计方差贡献 85% 的分量数 fprintf(保留 %d 个成分可解释 %.1f%% 方差\n, k, cumExp(k));coeff是载荷矩阵每列一个主成分score是投影后的得分直接可以拿去训练分类器latent是各成分方差explained是方差贡献百分比。NumComponents不指定时全部返回指定后explained仍然按全部成分给百分比这点在算累计贡献时要留意。VariableWeights,variance等价于先对每列做标准化再 PCA量纲单位不同的特征比如身高用厘米、体重用公斤一定要打开否则方差大的变量会独占第一主成分。因子分析和 PCA 的区别在于前者假设观测变量由少数潜在因子加噪声生成用factoran需要指定因子个数并选择旋转方式rotate可取varimax、promax。PCA 是无模型的方差分解因子分析有明确的测量模型解释潜在结构时后者更合适纯做压缩和去相关用 PCA 就够。方法目标关键输出什么时候用pca最大化投影方差载荷、得分、方差贡献去相关、降维、可视化factoran解释潜变量结构载荷、因子方差、特有方差量表分析、潜变量建模特征选择保留原始变量变量子集与重要性需要可解释系数5. 用 bootstrap 与并行把统计仿真跑到生产可用5.1 bootci 与 jackknife 的区间估计样本量不大、分布未知时基于正态近似的置信区间往往盖不住真值。bootci用重抽样直接构造经验区间比公式法稳健得多。rng(7); % 固定种子保证区间可复现 x exprnd(3, 500, 1); ci bootci(2000, mean, x, ... Alpha, 0.05, ... Type, per); % 百分位法 fprintf(均值 95%% 置信区间[%.3f, %.3f]\n, ci(1), ci(2));第一个参数是重抽样次数2000 次通常够用要求尾部精度时提到 10000。第二个参数是统计量句柄可以是mean、median、std也可以是自己写的匿名函数。Type取per是百分位法取bca是偏差校正加速法后者对小样本和偏态分布明显更准代价是计算量大几倍取student需要额外提供方差的估计实际用得少。jackknife是留一法适合估计偏差和标准误但它假设统计量光滑遇到中位数这类不光滑统计量会失效此时仍然用 bootstrap。5.2 向量化与 parfor 的收益边界把重抽样次数从 2000 提到 20000单线程可能要跑几分钟。加速有两条路向量化和并行。向量化是把每次重抽样写成矩阵操作用randi一次性生成所有索引避免循环。nBoot 20000; n numel(x); idx randi(n, n, nBoot); % n×nBoot 的索引矩阵内存换时间 bootStat mean(x(idx), 1); % 沿第 1 维求均值一次算完这段代码的内存开销是n × nBoot × 8字节n500、nBoot20000 时约 76 MB可以接受n 到几万时就要分块否则会触发内存不足。分块的做法是按列切成若干批每批用一次mean结果拼起来。当单次重抽样的代价很高比如每次都要重新拟合一个模型向量化就无能为力这时用parfor。p parpool(Processes, 4); rng(0, combRecursive); % 主种子并行下推荐组合递归生成器 bootStat zeros(nBoot, 1); parfor b 1:nBoot s randsample(n, n, true); % 有放回重抽样 bootStat(b) mean(x(s)); % 每次迭代独立无跨迭代依赖 end delete(p);parfor的收益取决于每次迭代的计算量与通信量的比值。像上面这种只求一个均值并行反而可能比向量化慢因为启动池和传输数据的开销占了大头。收益明显的是每次迭代里要跑fitglm、fitlm或者一次蒙特卡洛定价这类耗时操作。rng在主进程设置后并行池里每个 worker 会从主流派生子流worker 数量固定时结果可复现worker 数一变结果就会变这一点写进脚本注释里比写在文档里有用。本文还有配套的精品资源点击获取
MATLAB数理统计高级篇：分布对象、假设检验与回归建模实战
简介面向需要系统掌握MATLAB高级数理统计功能的科研与工程人员这套压缩包聚焦实战应用内容涉及多变量分析、假设检验、非参数检验、回归与拟合、时间序列分析、随机过程、生存分析、贝叶斯统计以及聚类判别等核心主题可帮助读者突破常规统计建模瓶颈。包内共2个文件含1段mp4视频演示与1个m脚本代码文件视频便于跟随讲解逐步理解统计方法的实现流程脚本可按自己的数据结构直接修改运行从而加深对代码细节的掌握压缩包整体约15.23MB轻量便于下载。目前已有115人学习下载适合希望在数据分析、预测建模与科学计算中提升统计应用能力的读者。文件结构简洁m脚本中的注释与视频讲解相互配合便于理解和复用通过这些内容可同时获得方法讲解、可执行代码与统计思维训练并能快速上手PCA、ARIMA、MCMC等实用技术。1. 从「11 matlab数理统计高级篇.zip」说起一个压缩包背后要补的三块硬功夫拿到一个写着「11 matlab数理统计高级篇.zip」的资源包多数人的第一反应是解压、翻目录、找看起来最像「高级」的那个 .m 文件然后直接 F5 运行。跑不通的占大多数报错信息往往还都很朴素变量未定义、矩阵维度不一致、函数或变量找不到。真正卡住人的从来不是代码本身而是压缩包里默认你已经会的那三块东西——概率分布对象怎么构造和复用、统计检验在不同数据结构下该选哪一个、仿真结果凭什么算可复现。看到 zip 先确认一件事它通常只交付脚本不交付随机种子的约定、不交付参数估计的置信区间口径、也不交付多重比较的校正方法。这个包要能用起来缺的正是这三块硬功夫而不是更多函数名。2. MATLAB 数理统计的分布对象体系与可复现随机数2.1 用 makedist 与 fitdist 把分布变成可传递的对象很多从入门教程过渡到 matlab 数理统计高级篇的人会一直停在normpdf、normcdf、normrnd这类「函数名 分布名前缀」的写法上。这套写法的致命问题是分布一旦换掉整套代码要逐行改前缀分布参数一旦从已知变成待估函数签名也对不上。较新的 MATLAB 版本提供的是面向对象的分布对象体系分布本身变成一个可以赋值、可以传参、可以当结构体字段存下来的变量。% 构造一个参数已知的分布对象后续所有函数复用同一对象 pd makedist(Normal, mu, 0, sigma, 1); % 从样本反推参数fitdist 返回一个已经拟合好的分布对象 x wblrnd(2, 1.5, 2000, 1); % 生成 Weibull 样本A2 B1.5 pdFit fitdist(x, Weibull); % 极大似然估计 disp(pdFit.ParameterNames); % 输出 {A,B}A 尺度B 形状 disp(pdFit.ParameterValues); % 输出估计值 ci paramci(pdFit, Alpha, 0.05); % 参数的 95% 置信区间 disp(ci); % 两列分别对应 A、B 的下上界这段代码做了三件事。makedist负责在参数已知时构造理论分布返回的对象可以直接喂给random、cdf、icdf、pdf也可以作为假设检验里CDF参数的对照。fitdist负责参数未知时的点估计默认走极大似然Distribution位置参数支持Normal、Weibull、Lognormal、Gamma、Kernel等一批名字。paramci是拟合后对象的方法Alpha默认 0.05Parameter还可指定只算某一个参数的区间。提示fitdist在样本被截断时一定要传Censoring否则极大似然会把删失点当成完整观测参数会系统性偏小。Kernel是里面唯一的非参数选项需要额外指定Kernel如epanechnikov、normal和Supportunbounded、positive、unit interval。它的ParameterValues是带宽和核中心不能拿去算 AIC这点在写脚本时最容易踩。2.2 rng 的状态控制与并行环境下的随机源选择数理统计里最容易被忽略、事后最难补救的一环是随机数能不能复现。同一份脚本在别人的机器上跑出不一样的 p 值几乎全部来自随机流没有固定。MATLAB 的rng提供三种典型用法对应三种不同目的。s rng(20240101, twister); % 固定种子并返回当前生成器状态结构体 a randn(1, 5); rng(s); % 把状态还原后续序列完全重放 b randn(1, 5); disp(isequal(a, b)); % 输出 1说明复现成功 rng(shuffle); % 需要每次不同时使用按时间做种子rng第一个参数是种子第二个参数是生成器类型常用取值有twisterMersenne Twister单线程默认、combRecursive组合递归周期极长适合并行、philox和threefry计数器型显式支持并行子流。保存返回的结构体再还原比只记一个种子更可靠因为状态结构体同时记录了生成器和位置。注意在parfor内部直接调用randnMATLAB 会给每个 worker 分配独立的子流以保证统计独立但这意味着结果会随 worker 数量变化。要严格复现先把 worker 数固定再在parfor之前用rng(seed,combRecursive)统一设置主种子。2.3 拟合优度诊断QQ 图、K-S 检验与 AIC 横向比选拿到一组样本后先问「它服从什么分布」这是 matlab 数理统计里出现频率最高的操作。凭直方图目测形状最容易出错尤其是尾部。工程上常用的判断顺序是先画 QQ 图看整体偏离再用 K-S 或 Anderson-Darling 做定量检验最后用 AIC 在若干候选分布之间比选。% 单样本 K-S检验样本是否来自指定理论分布 pd makedist(Normal, mu, mean(x), sigma, std(x)); [h, p, ksstat] kstest(x, CDF, pd); fprintf(K-S h%d p%.4g stat%.4f\n, h, p, ksstat); % 多候选分布的 AIC 比选 names {Normal, Weibull, Lognormal, Gamma}; aic zeros(numel(names), 1); for k 1:numel(names) pdK fitdist(x, names{k}); nll negloglik(pdK); % 负对数似然 aic(k) 2 * nll 2 * numel(pdK.ParameterValues); % AIC 2nll 2k end [~, best] min(aic); fprintf(最优分布%s\n, names{best});kstest的三个返回值分别是拒绝原假设的标志、p 值和检验统计量CDF传分布对象即可对照理论分布h1表示在默认 0.05 水平下拒绝。AIC 部分的关键是negloglik与参数个数的对应关系Normal和Lognormal各两个参数Weibull和Gamma也是两个所以这里自由度一致比的是拟合优度本身。候选分布典型适用场景需要留意的边界Normal测量误差、均值型指标对偏态和厚尾极敏感Weibull寿命、强度、失效时间形状参数小于 1 时尾部很重Lognormal收入、粒径、细胞数取对数后要再做一次正态检验Gamma等待时间、降水量、保险赔付尺度参数解释依赖业务背景提示K-S 检验对参数是用样本估计出来的情况偏保守p 值会偏大。样本量不大时优先用adtest或lillietest也可以直接做参数 bootstrap 的 K-S。3. 假设检验与方差分析从 ttest 到 anovan 的参数怎么设3.1 检验选型一张表定下 ttest2、ranksum 还是 signrank检验函数选错后面结论全部作废而且脚本不会报任何错。选型只看三件事比较的是均值还是分布位置、样本是否配对、方差是否齐性。把这些条件列成表比记函数名有用得多。场景函数调用必须设置的参数单样本均值对常数ttest(x, mu)Tail取both/right/left双样本、方差齐ttest2(x, y)Vartype,equal双样本、方差不齐ttest2(x, y)Vartype,unequalWelch 校正配对样本ttest(x, y)直接传两个等长向量即为配对不满足正态ranksum/signrankmethod取exact或approximate多组独立、非参数kruskalwallisdisplay取on方差不齐时硬用equal在两组样本量差异大的情况下第一类错误率会明显膨胀。稳健做法是先跑一次方差齐性检验再分支代码结构如下。load(carbig.mat, MPG, Origin); us MPG(strcmp(Origin, USA)); jp MPG(strcmp(Origin, Japan)); [~, vp] vartest2(us, jp); % 方差齐性检验只取 p 值 if vp 0.05 vtype unequal; % 方差不齐走 Welch else vtype equal; end [h, p, ci, stats] ttest2(us, jp, Vartype, vtype, Alpha, 0.01); fprintf(h%d t%.3f df%.1f p%.4g CI[%.3f, %.3f]\n, ... h, stats.tstat, stats.df, p, ci(1), ci(2));vartest2的第二个输出才是 p 值第一个是拒绝标志写成[h, vp]再取vp最不容易出错。ttest2返回的stats里有tstat和dfWelch 情形下df是小数这是正常现象。Alpha默认 0.05把它显式写出来是多组比较前控制族错误率的第一步。注意ttest(x, y)与ttest2(x, y)在语法上完全兼容传参错了也不会报错只会给出错误的结论。配对设计必须用ttest独立设计必须用ttest2。3.2 anovan 多因素方差分析与 multcompare 事后比较多因素设计的方差分析用anovan它比anova1、anova2笨一点因为要先把分组变量打包成 cell但换来的好处是支持任意因子数、不平衡设计和指定交互项。% y 为响应向量g1/g2/g3 为分组向量长度必须一致 [p, tbl, stats] anovan(y, {g1, g2, g3}, ... model, interaction, ... % 主效应 所有两两交互 varnames, {温度, 压力, 批次}, ... sstype, 3, ... % 不平衡数据用 III 型平方和 display, on); % 事后两两比较Bonferroni 控制族错误率 [c, m, h, gnames] multcompare(stats, ... CType, bonferroni, ... Alpha, 0.05, ... Display, on);model是这里最关键的参数。linear只做主效应interaction加所有两两交互full加全部交互项也可以直接传一个矩阵来自定义哪些因子参与交互比如[1 0 0; 0 1 0; 0 0 1; 1 1 0]表示三个主效应加上「因子 1 与因子 2」的交互。sstype默认是 3样本量在各组间不均衡时I 型和 II 型平方和会随因子进入顺序变化III 型不会所以不平衡设计务必显式写 3。multcompare的CType决定事后比较的校正方式bonferroni最保守tukey-kramer均衡设计下更紧scheffe适合任意线性组合的比较lsd不做校正只适合事先已经锁定了少数几对比较的情况。返回的c矩阵前两列是要比较的两组编号末列是 p 值m是各组的估计边际均值。提示anovan要求分组变量要么是数值向量要么是 cell 字符串数组。如果原始数据是 table先用categorical转换再取double否则因子水平顺序会按字典序而不是你期望的业务顺序排。3.3 p 值多重校正手写 Benjamini-Hochberg 控制 FDR一次跑几十上百个检验不做校正就报「显著」是统计脚本里最严重的错误。Bonferroni 控制的是族错误率比较次数一多就几乎什么都检不出来工程上更常用的是 Benjamini-Hochberg 控制错误发现率。MATLAB 的统计与机器学习工具箱没有直接提供 BH 函数手写反而更透明。p [0.001, 0.008, 0.039, 0.041, 0.12, 0.33, 0.65]; m numel(p); [sorted_p, idx] sort(p(:)); % 升序排列并记录原位置 bh sorted_p .* m ./ (1:m); % 第 i 个乘 m/i bh flipud(cummin(flipud(bh))); % 从大到小取累计最小保证单调 q zeros(m, 1); q(idx) min(bh, 1); % 按原顺序还原并截断到 1 disp(q);三个步骤对应 BH 的定义。先把 p 值升序排列并记下原始下标第 i 小的 p 值乘以 m/i 得到校正后的 q 值然后从最大的一头往回取累计最小值这一步是为了让 q 值单调不减否则会出现「更小的原始 p 值算出更大的 q 值」这种反直觉结果。最后按原索引还原并截断到 1。判定显著时把 q 与目标 FDR 水平常取 0.05 或 0.10比而不是跟原始 p 值比。如果是基因表达这类场景Bioinformatics Toolbox 里的mafdr直接支持BHFDR选项还额外实现了 Storey 的 q 值方法估计 π0 后检验效力更高。两者结论方向一致但 q 值数值会有差异写报告时要注明用的是哪一种。4. 回归与降维fitlm、fitglm、pca 的工程化用法4.1 fitlm 稳健回归与残差诊断量线性模型在 matlab 数理统计高级篇里出现的频次仅次于假设检验但多数脚本只打印mdl.Coefficients就收工不看残差也不看杠杆点。fitlm的正确用法是把 table 直接传进去让变量名自动进入公式然后用诊断图确认假设是否成立。tbl table(x1, x2, x3, y, VariableNames, {x1,x2,x3,y}); mdl fitlm(tbl, y ~ x1 x2 x3 x1:x2, ... RobustOpts, on, ... % bisquare 稳健拟合抑制离群点 Intercept, true); disp(mdl.Coefficients); % Estimate / SE / tStat / pValue disp(mdl.Rsquared.Adjusted); % 调整 R²比 R² 更可信 anova(mdl); % 各项平方和与 F 检验 plotResiduals(mdl, fitted); % 残差对拟合值看是否喇叭形 plotDiagnostics(mdl, leverage); % 杠杆值找强影响点RobustOpts取on时用 bisquare 权重迭代重加权等价于对离群点降权如果想自定义权重函数传结构体struct(wfun,fair,tune,1.5)即可tune是调优常数。Rsquared是结构体含Ordinary、Adjusted、Ordinary三个字段多元回归里看Adjusted。残差图出现明显的漏斗形说明异方差此时应改用稳健标准误或对响应做变换而不是直接读 p 值。逐步回归用stepwiselmPEnter和PRemove默认分别是 0.05 和 0.10后一个比前一个大是为了防止变量在边界上反复进出。这个函数对样本量小、候选变量多的数据容易过拟合用它筛出来的模型一定要在留出集上重新验证一次。4.2 fitglm 做二分类与链接函数选择响应是 0/1 或者计数时不能再用fitlm改用fitglm。它比旧的glmfit友好在于直接接受 table 和公式并且自带predict和devianceTest。% X 是 n×p 预测矩阵y 是 n×1 的 0/1 向量 mdl fitglm(X, y, Distribution, binomial, ... Link, logit, ... % 也可取 probit、comploglog LikelihoodPenalty, jeffreys-prior); % 完全分离时稳定估计 prob predict(mdl, Xnew); % 输出属于正类的概率 label prob 0.5; % 模型整体显著性与常数模型做偏差检验 devianceTest(mdl);Distribution取binomial时Link可选logit、probit、comploglog取poisson时通常配log。LikelihoodPenalty是较新版本引入的正则化选项jeffreys-prior在样本量小、类别完全可分的情况下能显著稳定系数否则标准误会膨胀到不可用。predict返回的概率不要直接用 0.5 切类别不平衡时应该按业务代价调整阈值或者用 ROC 曲线找约登指数最大的点。如果坚持用glmfit要自己给设计矩阵补一列全 1返回的b是系数、dev是偏差、stats.p是每个系数的 p 值接口更底层但输出更杂。新脚本建议直接用fitglm。4.3 pca 降维与因子分析在特征压缩中的落地维度高了以后回归系数不稳定先降维再建模是常见流程。pca的输出字段多写脚本时按名字取不要依赖位置顺序。[coeff, score, latent, tsquared, explained, mu] pca(X, ... NumComponents, 10, ... % 只保留前 10 个主成分 Algorithm, svd, ... % svd 精度最高eig 更快 Centered, true, ... % 默认按列中心化 VariableWeights, variance); % 变量量纲差异大时按方差加权 cumExp cumsum(explained); k find(cumExp 85, 1); % 取到累计方差贡献 85% 的分量数 fprintf(保留 %d 个成分可解释 %.1f%% 方差\n, k, cumExp(k));coeff是载荷矩阵每列一个主成分score是投影后的得分直接可以拿去训练分类器latent是各成分方差explained是方差贡献百分比。NumComponents不指定时全部返回指定后explained仍然按全部成分给百分比这点在算累计贡献时要留意。VariableWeights,variance等价于先对每列做标准化再 PCA量纲单位不同的特征比如身高用厘米、体重用公斤一定要打开否则方差大的变量会独占第一主成分。因子分析和 PCA 的区别在于前者假设观测变量由少数潜在因子加噪声生成用factoran需要指定因子个数并选择旋转方式rotate可取varimax、promax。PCA 是无模型的方差分解因子分析有明确的测量模型解释潜在结构时后者更合适纯做压缩和去相关用 PCA 就够。方法目标关键输出什么时候用pca最大化投影方差载荷、得分、方差贡献去相关、降维、可视化factoran解释潜变量结构载荷、因子方差、特有方差量表分析、潜变量建模特征选择保留原始变量变量子集与重要性需要可解释系数5. 用 bootstrap 与并行把统计仿真跑到生产可用5.1 bootci 与 jackknife 的区间估计样本量不大、分布未知时基于正态近似的置信区间往往盖不住真值。bootci用重抽样直接构造经验区间比公式法稳健得多。rng(7); % 固定种子保证区间可复现 x exprnd(3, 500, 1); ci bootci(2000, mean, x, ... Alpha, 0.05, ... Type, per); % 百分位法 fprintf(均值 95%% 置信区间[%.3f, %.3f]\n, ci(1), ci(2));第一个参数是重抽样次数2000 次通常够用要求尾部精度时提到 10000。第二个参数是统计量句柄可以是mean、median、std也可以是自己写的匿名函数。Type取per是百分位法取bca是偏差校正加速法后者对小样本和偏态分布明显更准代价是计算量大几倍取student需要额外提供方差的估计实际用得少。jackknife是留一法适合估计偏差和标准误但它假设统计量光滑遇到中位数这类不光滑统计量会失效此时仍然用 bootstrap。5.2 向量化与 parfor 的收益边界把重抽样次数从 2000 提到 20000单线程可能要跑几分钟。加速有两条路向量化和并行。向量化是把每次重抽样写成矩阵操作用randi一次性生成所有索引避免循环。nBoot 20000; n numel(x); idx randi(n, n, nBoot); % n×nBoot 的索引矩阵内存换时间 bootStat mean(x(idx), 1); % 沿第 1 维求均值一次算完这段代码的内存开销是n × nBoot × 8字节n500、nBoot20000 时约 76 MB可以接受n 到几万时就要分块否则会触发内存不足。分块的做法是按列切成若干批每批用一次mean结果拼起来。当单次重抽样的代价很高比如每次都要重新拟合一个模型向量化就无能为力这时用parfor。p parpool(Processes, 4); rng(0, combRecursive); % 主种子并行下推荐组合递归生成器 bootStat zeros(nBoot, 1); parfor b 1:nBoot s randsample(n, n, true); % 有放回重抽样 bootStat(b) mean(x(s)); % 每次迭代独立无跨迭代依赖 end delete(p);parfor的收益取决于每次迭代的计算量与通信量的比值。像上面这种只求一个均值并行反而可能比向量化慢因为启动池和传输数据的开销占了大头。收益明显的是每次迭代里要跑fitglm、fitlm或者一次蒙特卡洛定价这类耗时操作。rng在主进程设置后并行池里每个 worker 会从主流派生子流worker 数量固定时结果可复现worker 数一变结果就会变这一点写进脚本注释里比写在文档里有用。本文还有配套的精品资源点击获取
