简介：一份基于研究生教材《数理统计》例4.4.1编写的多元线性回归及显著性检验Matlab程序文档，面向统计学习者、数据分析人员和需要快速完成回归建模的开发者。文档完整呈现了程序原理、数据存储格式与可直接运行的Matlab代码，并在教材原有回归方程显著性分析基础上，单独对每个自变量的回归系数进行t检验，同时支持用户自定义显著性水平α，自动剔除不显著变量，提高计算精度。程序采用最小二乘法β=(X'X)^(-1)X'Y估计回归系数，依次完成总平方和、残差平方和计算，并利用F统计量判断回归方程整体显著性，再通过t统计量与临界值比较评估各系数是否显著不为零。资源包仅含1个docx文件，大小约50KB，无需解压即可在Word中查看源码与说明。已有246人学习下载，适合需要借助Matlab理解最小二乘估计、F检验与逐步回归筛选流程的读者；替换Excel数据后还可应用于其他多元数据集，是教学演示与科研实践的双重便利工具。
1. 为什么教材例题要改成交互式显著性检验程序
多元线性回归的显著性检验，在很多教材里只讲一次F检验和一轮t检验，剔除一个变量后就戛然而止。研究生教材《数理统计》例4.4.1就是这样，删除x1之后没有再对x2、x3做回归系数显著性检验。但实际回归诊断要求“剔除一个变量后重新拟合、重新检验”，因为删变量会改变剩余系数的估计值和标准误，原来的t值全部失效。这个Matlab程序把循环补齐了，并在每次迭代中输出cii、临界边界和系数，让使用者看清每一步变量筛选的依据。程序还允许用户输入任意显著性水平α，而不是写死0.05。对互联网行业的数据分析岗位来说，做特征筛选和归因分析时经常要回答“这个变量到底有没有解释力”，这类逐步检验程序比直接套Excel回归工具更可控，也更容易嵌入论文或实验报告。
2. 最小二乘估计与回归方程显著性F检验的Matlab实现
2.1 从Excel到设计矩阵X
先看数据读取和设计矩阵构造：
data = xlsread('jc_p133_example.xls', 'sheet1'); xi = data(:, 1:end-1); [n, k] = size(data); k = k - 1; X = [ones(n, 1) xi]; Y = data(:, end);xlsread直接读取Excel数值区域，要求工作表中不能有表头文本。如果带表头，需要改成xlsread('...', 'sheet1', 'A2:F50')，否则读进来是NaN。xi取前k列自变量，Y取最后一列因变量。size(data)返回[n, k+1]，所以真正自变量个数是k-1，这里n是样本量。X = [ones(n,1) xi]是在自变量矩阵前插一列全1，对应截距项β0。
最小二乘求解用的是左除：
beta_mao = ((X' * X) \ X' * Y)';左除\在Matlab内部会做QR分解或Cholesky分解，数值稳定性比inv(X'*X)*X'*Y好得多。当自变量存在一定相关性时，直接用inv可能得到很大的对角元素，而左除能抑制部分浮点误差。程序把结果转置成行向量，后面打印β时方便按顺序输出。
程序里的核心变量可以先用一张表理清：
| 变量 | 维度 | 含义 |
|---|---|---|
data | n×(k+1) | 原始数值矩阵，最后一列是y |
xi | n×k | 自变量矩阵，不含截距 |
X | n×(k+1) | 补全1列后的设计矩阵 |
beta_mao | 1×(k+1) | 最小二乘回归系数 |
St_square | 标量 | 总平方和SST |
Sr_square | 标量 | 回归平方和SSR |
Se_square | 标量 | 残差平方和SSE |
cii | 1×(k+1) | (X'X)^-1对角线元素 |
ci | 1×k | 系数显著性临界边界 |
2.2 平方和分解与回归平方和的计算
接下来是平方和计算：
x_ba = mean(xi); y_ba = mean(Y); St_square = sum(Y.^2) - n * y_ba^2; lxy = sum((xi - ones(n,1)*x_ba) .* ((Y - y_ba) * ones(1,k))); Sr_square = sum(beta_mao(2:end) .* lxy); Se_square = St_square - Sr_square;St_square是总平方和SST，等于Σy_i² - n·ȳ²，也就是Σ(y_i-ȳ)²。这里先展开平方和再减均值项，能少一次矩阵减法。lxy是自变量的中心化矩阵与因变量中心化向量的交叉积，结果是一个1×k的行向量，第j个元素是Σ(x_ij-x̄_j)(y_i-ȳ)。回归平方和SSR直接用系数与lxy的内积得到：SSR = Σβ_j·lxy_j。这么做的好处是不需要先计算拟合值ŷ，也避免了构造n×n投影矩阵，样本量较大时内存优势明显。
注意beta_mao(2:end)剔除了截距项，因为截距不参与回归平方和。Se_square通过减法得到，理论上SSE=SST-SSR。如果数据严重多重共线性，SSR可能略大于SST，需要检查原始数据。我一般会在这里加一句判断：
assert(Se_square > 0, '残差平方和为负，可能存在严重多重共线性');2.3 F检验与临界值折算
回归方程显著性检验的代码：
F_fenweidian = finv(1 - F_alpha, k, n - k - 1); c = k / (n - k - 1) * F_fenweidian; if Sr_square / Se_square > c fprintf('拒绝H0，回归方程显著\n'); else fprintf('接受H0，回归方程不显著\n'); end这段代码初看容易懵：为什么不是标准的(Sr/k)/(Se/(n-k-1))？程序直接比较Sr/Se和一个调整后的临界值c，实际上是对不等式做了等价变形：
(Sr/k)/(Se/(n-k-1)) > F_{1-α}(k, n-k-1)
等价于Sr/Se > (k/(n-k-1))·F_{1-α}(k, n-k-1)。
所以c是折算后的F临界值，Sr/Se是简化统计量。finv(1-F_alpha, k, n-k-1)取得的是下侧分位数，由于F检验是右侧单尾，下侧概率要取1-α。例如α=0.02，k=3，n=14时，finv(0.98,3,10)约等于多少？程序会计算。如果把1-F_alpha误写成F_alpha，临界值会变小很多，导致错误拒绝H0。这是使用finv最常见的坑。
3. 回归系数显著性检验的循环剔除机制与实现细节
3.1 用F分布临界值代替t检验
回归系数β_j的显著性通常用t检验，统计量t_j = β_j / sqrt(cii(j)·Se/(n-k-1))。因为t_j²服从F(1, n-k-1)，程序直接用F分布的分位数构造系数临界值ci：
F_fenweidian_1 = finv(1 - F_alpha, 1, n - k - 1); ci = sqrt(cii(2:end) * Se_square * F_fenweidian_1 / (n - k - 1));cii是inv(X'*X)的对角线，cii(2:end)去掉截距对应的第一项。ci的每个元素表示“当前α下该系数绝对值的临界边界”。如果某个|β_j| > ci(j)，说明该系数显著不为零；否则不显著，应当从模型中剔除。
这里有一个容易误解的地方：ci并不是标准误，而是标准误乘以sqrt(F_fenweidian_1)。因为t临界值t_{1-α/2}的平方等于F_{1-α}(1, df)，所以ci相当于“最小显著系数值”。
3.2 定位最不显著变量
程序在每轮循环中要找最不显著的那个变量：
fi_xin = beta_1tok.^2 ./ cii(1:end-1)'; min_fi = min(fi_xin); beta_index = find(fi_xin == min_fi) + 1;fi_xin保存的是每个变量对应的t²值（等价于F值），取最小值就找到了最不显著的变量。find(fi_xin == min_fi)在多个变量数值完全相同时返回多个索引，此时程序取第一个。这在真实数据里不常见，但如果是人造数据，可能同时有多个完全相等的t²，建议用find(...,1)只取第一个。
这里要注意cii(1:end-1)和beta_1tok的长度都等于当前保留的变量个数，而beta_index是相对beta_mao的索引，所以要+1。因为beta_mao第一位是截距β0，后面才是β1、β2、β3。
3.3 删除变量后的系数快速更新
找到最不显著变量并标记为删除后，程序做了一次巧妙的系数更新：
beta_mao = beta_mao - beta_mao(beta_index) / cii(beta_index) * cij(beta_index, :); beta_mao(beta_index) = [];这段代码等价于重新用新设计矩阵做一次最小二乘，但避免了重新求逆。推导思路：设A = X'X，b = X'y，当前解β = A^{-1}b。删除第j个变量相当于把A的第j行第j列去掉构成子矩阵A_j，对应的解β_j = A_j^{-1} b_j。利用分块矩阵求逆公式，可以得到更新式：新的系数向量等于旧系数减去β_j / A_jj乘以A^{-1}的第j行（除第j个元素外）。程序里cij(beta_index,:)就是A^{-1}的第beta_index行，cii(beta_index)是对角元。更新后把该位置系数删除，再同步删除X中对应列。
这个技巧在变量上百个时的加速效果明显。但要注意，更新后Se_square也会变，必须用新X重新计算。程序在循环末尾重新计算Sr_square和Se_square，这部分不能省。
3.4 可读性输出与索引保持
程序输出回归系数用了动态格式字符串：
fmt_str1 = ''; for i1 = 2:k+1 fmt_str1 = [fmt_str1 'β' num2str(i1-1+num_of_loop) ' = %0.4f\r']; end fprintf(['β0 = %0.4f\r' fmt_str1], beta_mao);这里的编号i1-1+num_of_loop是为了让输出对应原始变量编号。因为每次删除变量后，beta_mao里的位置会缩短，但程序中变量原本是x1、x2、x3，所以用num_of_loop补偿已经删除的个数。直接看这段代码会有点绕，我一般会用一个单独的kept_index向量记录当前剩余变量的原始编号，输出时直接用原始编号，避免加法混乱。
4. 数据组织、α参数与程序移植性设计
4.1 Excel数据格式与预处理
程序对Excel数据的要求是“按x1,x2,…,xk,y列存放”，不带表头。原始文本粘贴时，数据之间有多个空格，Excel的“数据-分列-按空格”可以转成列。注意第一列第一个数如果是变量值，没问题。下面是一个简单的示例格式：
| x1 | x2 | x3 | y |
|---|---|---|---|
| ... | ... | ... | ... |
如果原始数据里出现空行或空单元格，xlsread默认填NaN，回归前必须检查。我会在读取后加：
if any(isnan(data(:))) error('数据中包含NaN，请检查Excel表格'); end否则X'*X里出现NaN，后面所有检验都会静默失败或给出错误结论。
4.2 显著性水平α的交互输入与自动化改造
程序用input交互获取α，并做了输入校验：
F_alpha = input('请输入显著性水平α(0<α<1): '); while ~(isscalar(F_alpha) && F_alpha < 1 && F_alpha > 0) F_alpha = input('输入有误，请重新输入α: '); end在脚本里这样没问题，但如果你要在循环里跑100个α，input会卡住。我通常改成：
if ~exist('F_alpha', 'var') F_alpha = 0.05; % 默认值 end或者把整个程序封装成function reg_analysis(data, alpha)，把F_alpha作为参数传入。这样既能交互，又能批量跑：
alpha_list = [0.01 0.02 0.05 0.10]; for a = alpha_list reg_analysis(data, a); end4.3 移植到其他数据集时的三个检查点
第一，样本量必须大于自变量个数加1，即n > k+1，否则X'X奇异，左除会给出警告，cii可能为负数。第二，不能有完全共线性列。比如某一列刚好是另一列的2倍，X'X行列式为0，系数无法唯一估计。用rcond(X'*X)可以快速判断：
if rcond(X'*X) < 1e-12 warning('设计矩阵接近奇异，请检查自变量相关性'); end第三，数据量纲不要差太大。比如一列是几十，另一列是几万，X'*X的条件数会很大，虽然左除比inv稳定，但cii数值仍可能不稳定。常见的做法是标准正态化后再回归，标准化后的系数是标准化系数，解释时要回退。程序没有做标准化，所以输入原始数据时最好先做一次无量纲化。
4.4 与Matlab内置工具的差异
Matlab自带regress和fitlm也能做回归和显著性检验，但它们是“一次性”检验，不会自动逐步剔除不显著变量。stepwiselm可以做逐步回归，但它的判据是AIC或BIC，不是自定义α下的F检验。这个程序的价值在于完全透明：每一次剔除都打印cii和临界值，用户可以核对每个变量的t²到底是多少。这对教学和分析报告的附录很有用，内置函数做不到这种中间状态输出。
5. 用α=0.01和α=0.02的输出验证检验流程与边界
5.1 两次运行的差异对照
同一份数据，α=0.01时程序依次剔除了x1、x2，最后保留x3；α=0.02时三个变量都保留。差异直观对比如下：
| α | 回归方程检验 | 系数检验结果 |
|---|---|---|
| 0.01 | 拒绝H0，方程显著 | 剔除x1后剔除x2，最终保留x3 |
| 0.02 | 拒绝H0，方程显著 | x1、x2、x3全部显著 |
这说明显著性水平越严格，剔除的变量越多。α=0.01意味着要求证据更强，p值必须小于0.01才能保留变量。α=0.02虽然不是常用的0.05，但在比较试验中可以看到阈值变化对结果的影响。实际业务中，我一般先跑α=0.05，再跑α=0.01，如果结论差异大，就说明数据中存在边界显著的变量，需要结合业务判断要不要保留。
5.2 用regress函数交叉验证
程序输出的结论可以通过内置函数快速验证：
[b, ~, ~, ~, stats] = regress(Y, X); fprintf('R²=%.4f, F=%.4f, p=%.6f\n', stats(1), stats(2), stats(3));stats向量的第2、3个分别是F统计量和p值。如果p值小于α，则“拒绝H0，方程显著”的结论与程序一致。对于系数显著性，可以用fitlm查看每个系数的p值：
mdl = fitlm(xi, Y); disp(mdl.Coefficients);mdl.Coefficients里有每个变量的tStat和pValue，和程序计算的fi_xin对应的p值应当一致。不一致时先检查cii是否计算正确，再看自由度是否用了n-k-1而不是n-k。
5.3 α作为外部参数批量扫描的技巧
最后分享一个很实用的技巧。把程序主体抽成一个函数后，可以用arrayfun或for循环批量扫描α：
alpha_list = 0.01:0.005:0.05; for idx = 1:length(alpha_list) fprintf('\n======== α = %.3f ========\n', alpha_list(idx)); reg_analysis(data, alpha_list(idx)); end这样一次运行就能看到不同显著性水平下变量筛选的变化，很适合做敏感性分析。输出量大时用diary记录：
diary('alpha_sensitivity.txt'); diary on; % 批量跑 diary off;diary会把所有fprintf输出重定向到文件，论文或汇报里可以直接引用。注意diary记录的是代码运行时的完整输出，包括之前的历史输出，所以在diary on之前先clear命令窗口比较干净。
本文还有配套的精品资源，点击获取
