简介一份基于研究生教材《数理统计》例4.4.1编写的多元线性回归及显著性检验Matlab程序文档面向统计学习者、数据分析人员和需要快速完成回归建模的开发者。文档完整呈现了程序原理、数据存储格式与可直接运行的Matlab代码并在教材原有回归方程显著性分析基础上单独对每个自变量的回归系数进行t检验同时支持用户自定义显著性水平α自动剔除不显著变量提高计算精度。程序采用最小二乘法β(XX)^(-1)XY估计回归系数依次完成总平方和、残差平方和计算并利用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, k1]所以真正自变量个数是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可能得到很大的对角元素而左除能抑制部分浮点误差。程序把结果转置成行向量后面打印β时方便按顺序输出。程序里的核心变量可以先用一张表理清变量维度含义datan×(k1)原始数值矩阵最后一列是yxin×k自变量矩阵不含截距Xn×(k1)补全1列后的设计矩阵beta_mao1×(k1)最小二乘回归系数St_square标量总平方和SSTSr_square标量回归平方和SSRSe_square标量残差平方和SSEcii1×(k1)(XX)^-1对角线元素ci1×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通过减法得到理论上SSESST-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.02k3n14时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分布的分位数构造系数临界值ciF_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 XXb Xy当前解β 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:k1 fmt_str1 [fmt_str1 β num2str(i1-1num_of_loop) %0.4f\r]; end fprintf([β0 %0.4f\r fmt_str1], beta_mao);这里的编号i1-1num_of_loop是为了让输出对应原始变量编号。因为每次删除变量后beta_mao里的位置会缩短但程序中变量原本是x1、x2、x3所以用num_of_loop补偿已经删除的个数。直接看这段代码会有点绕我一般会用一个单独的kept_index向量记录当前剩余变量的原始编号输出时直接用原始编号避免加法混乱。4. 数据组织、α参数与程序移植性设计4.1 Excel数据格式与预处理程序对Excel数据的要求是“按x1,x2,…,xk,y列存放”不带表头。原始文本粘贴时数据之间有多个空格Excel的“数据-分列-按空格”可以转成列。注意第一列第一个数如果是变量值没问题。下面是一个简单的示例格式x1x2x3y............如果原始数据里出现空行或空单元格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 k1否则XX奇异左除会给出警告cii可能为负数。第二不能有完全共线性列。比如某一列刚好是另一列的2倍XX行列式为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最终保留x30.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命令窗口比较干净。本文还有配套的精品资源点击获取 基于Matlab的多元线性回归逐步显著性检验程序实现
简介一份基于研究生教材《数理统计》例4.4.1编写的多元线性回归及显著性检验Matlab程序文档面向统计学习者、数据分析人员和需要快速完成回归建模的开发者。文档完整呈现了程序原理、数据存储格式与可直接运行的Matlab代码并在教材原有回归方程显著性分析基础上单独对每个自变量的回归系数进行t检验同时支持用户自定义显著性水平α自动剔除不显著变量提高计算精度。程序采用最小二乘法β(XX)^(-1)XY估计回归系数依次完成总平方和、残差平方和计算并利用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, k1]所以真正自变量个数是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可能得到很大的对角元素而左除能抑制部分浮点误差。程序把结果转置成行向量后面打印β时方便按顺序输出。程序里的核心变量可以先用一张表理清变量维度含义datan×(k1)原始数值矩阵最后一列是yxin×k自变量矩阵不含截距Xn×(k1)补全1列后的设计矩阵beta_mao1×(k1)最小二乘回归系数St_square标量总平方和SSTSr_square标量回归平方和SSRSe_square标量残差平方和SSEcii1×(k1)(XX)^-1对角线元素ci1×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通过减法得到理论上SSESST-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.02k3n14时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分布的分位数构造系数临界值ciF_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 XXb Xy当前解β 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:k1 fmt_str1 [fmt_str1 β num2str(i1-1num_of_loop) %0.4f\r]; end fprintf([β0 %0.4f\r fmt_str1], beta_mao);这里的编号i1-1num_of_loop是为了让输出对应原始变量编号。因为每次删除变量后beta_mao里的位置会缩短但程序中变量原本是x1、x2、x3所以用num_of_loop补偿已经删除的个数。直接看这段代码会有点绕我一般会用一个单独的kept_index向量记录当前剩余变量的原始编号输出时直接用原始编号避免加法混乱。4. 数据组织、α参数与程序移植性设计4.1 Excel数据格式与预处理程序对Excel数据的要求是“按x1,x2,…,xk,y列存放”不带表头。原始文本粘贴时数据之间有多个空格Excel的“数据-分列-按空格”可以转成列。注意第一列第一个数如果是变量值没问题。下面是一个简单的示例格式x1x2x3y............如果原始数据里出现空行或空单元格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 k1否则XX奇异左除会给出警告cii可能为负数。第二不能有完全共线性列。比如某一列刚好是另一列的2倍XX行列式为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最终保留x30.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命令窗口比较干净。本文还有配套的精品资源点击获取