简介：面向需要系统掌握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 样本，A=2 B=1.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.05，'Parameter'还可指定只算某一个参数的区间。
提示：
fitdist在样本被截断时一定要传'Censoring'，否则极大似然会把删失点当成完整观测，参数会系统性偏小。
'Kernel'是里面唯一的非参数选项，需要额外指定'Kernel'（如'epanechnikov'、'normal'）和'Support'（'unbounded'、'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第一个参数是种子，第二个参数是生成器类型，常用取值有'twister'（Mersenne Twister，单线程默认）、'combRecursive'（组合递归，周期极长，适合并行）、'philox'和'threefry'（计数器型，显式支持并行子流）。保存返回的结构体再还原，比只记一个种子更可靠，因为状态结构体同时记录了生成器和位置。
注意：在
parfor内部直接调用randn，MATLAB 会给每个 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'传分布对象即可对照理论分布；h=1表示在默认 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','unequal'（Welch 校正） |
| 配对样本 | ttest(x, y) | 直接传两个等长向量即为配对 |
| 不满足正态 | ranksum/signrank | 'method'取'exact'或'approximate' |
| 多组独立、非参数 | kruskalwallis | 'display'取'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和df，Welch 情形下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 值。
逐步回归用stepwiselm，'PEnter'和'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字节，n=500、nBoot=20000 时约 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 数一变结果就会变，这一点写进脚本注释里比写在文档里有用。
本文还有配套的精品资源，点击获取
