简介面向2023年美国大学生数学建模竞赛C题第一问提供一套基于Matlab的完整解决方案特别适合正在备赛的建模选手以及需要处理类似文本数据的科研人员。压缩包共五个文件包括两个xlsx表格分别存放原始词频数据与处理结果两个m脚本承担初始化与主流程求解任务一个mat数据文件保存了中间变量整体大小约67KB轻量易读。代码覆盖从数据读取、清洗到模型构建、结果验证与可视化的完整链路并针对Wordle数据集做了日期与词频维度的特征分析便于快速复现赛题图表。源码中附有必要的注释模块划分清晰学习者可对照赛题逐步调试。目前已有1008人学习下载说明该方案具有较高的参考价值。读者通过研读源码既能掌握Matlab处理真实题目的编程技巧也能迁移其中的数据分组统计思路至其他竞赛或科研场景中。1. 从原始数据到第一问的模型先明确这一问在解释什么2023 美赛 C 题给的是《纽约时报》Wordle 每天的玩家结果汇总列里是尝试 1 到 6 次的人数、废纸篓没猜出来的人数和总人数。第一问表面是“解释每日报告结果的变化”翻译成建模语言就是用一个数据驱动的方程把总人数随时间变化的趋势、按星期重复出现的波动以及某几天突然冒出来的异常值都解释清楚。核心不是把数字拟合得越准越好而是让模型能拆出可解释的成分。常见做法是先把数据读进来做探索性分析再用带周季节项的回归模型做主输出。MATLAB 很合适readtable 一行读 CSVdatetime 处理日期fitlm 求解很快。下面直接按能跑通的顺序写一遍源码运行环境是 MATLAB R2023b 及以上核心回归只依赖 Statistics Toolbox可选地使用 Econometrics Toolbox 画残差自相关Deep Learning Toolbox 用不到。适合准备美赛或国赛 C 题的参赛队以及想把时间序列回归做成可解释模型的工程师。2. 预处理与探索读入数据先看曲线长什么样2.1 用 readtable 和 datetime 把 Wordle 原始 CSV 读成可分析表拿到 2023 年 C 题的数据文件夹后最常见的问题是 Excel 打开看到的是乱码日期因为 CSV 里的日期是 M/d/yyyy 格式。用 MATLAB 读的时候不要用 importdata建议直接 readtable 然后自己控制日期格式clear; rng(42); % 读取原始数据文件路径改成你本机实际路径 T readtable(Wordle_Data_March_2022_to_Jan_2023.csv, ... TextType, string); % 日期列为 M/d/yyyy 格式要显式指定否则整列是 NaT T.Date datetime(T.Date, InputFormat, M/d/yyyy); T sortrows(T, Date); % 取当天总人数列原数据里这一列名通常叫 Total total T.Total_Number_of_Reported_Results; date T.Date; fprintf(样本范围: %s 到 %s, 共 %d 天\n, ... datestr(min(date)), datestr(max(date)), height(T));代码里先用TextType让字符串列以 string 类型进入内存避免后面字符数组拼接出问题datetime的InputFormat要和 CSV 里实际格式严格一致不然会得到一整列 NaT。很多参赛队在这里翻车Excel 把日期显示成 3/1/2022MATLAB 默认格式是 yyyy-MM-dd不加 InputFormat 会全部解析失败。读进来之后建议顺手做两层检查sum(isnat(date))必须为 0sum(isnan(total))必须为 0。这两个数只要有一个非零后面回归结果就不该继续往下走得先回去看是原始文件空行还是合并单元格引起的。2.2 用 7 天滑动平均将趋势和周内波动分离开原始总人数序列是一个高频时间序列直接看曲线只能看到“2022 年 4 月高、之后缓慢下降、偶尔有尖峰”看不出周末效应。我一般会先画一张原始曲线和 7 天滑动平均叠在一起的图% movmean 用 7 天窗口等价于消除周内周期 trend7 movmean(total, 7, Endpoints, discard); figure(Color, w, Position, [100 100 900 420]); plot(date, total, Color, [0.75 0.75 0.75]); hold on; plot(date(4:end-3), trend7, LineWidth, 2, ... Color, [0.85 0.33 0.10]); xlabel(日期); ylabel(当天报告总人数); legend({原始总人数, 7天滑动平均}, Location, best); grid on; box on;movmean的Endpoints指定两端点如何处理discard会让输出首尾各少 3 个观测。这样得到的 trend7 去掉的是以 7 天为周期的波动剩下的变化就主要是长期趋势加突发事件。对 2023 年 C 题的数据你会观察到日期区间原始现象建模时对应处理2022 年 3 月到 4 月中旬总人数快速上升分段上升趋势2022 年 4 月到 2023 年 1 月总人数稳步下滑分段下滑趋势2023 年 1 月初几天明显低点干预哑变量某些星期一或星期二出现低点周系数取负值这个表可以作为预处理的第一个输出放进论文评委看到的不只是一张图还能看出你已经把“趋势 周效应 异常”三个成分分开。2.3 对偏差过大的日期做 z-score 检查并记录干预点滑动平均能看趋势但没法客观地标出哪些日期算异常。更稳妥的做法是先算残差原始值减 trend7再用滚动标准差做 z-score。具体实现% 先用 NaN 把 movmean 丢弃的两端补回来再做逐点残差 resid total - [NaN(3,1); trend7; NaN(3,1)]; sd14 movstd(total, 14, Endpoints, discard); z abs(resid(4:end-3) ./ sd14); % 找出 z 超过 2.5 的日期 outlier_idx find(z 2.5); outlier_dates date(4:end-3); disp(table(outlier_dates(outlier_idx), total(4:end-3), z(outlier_idx), ... VariableNames, {日期, 总人数, z值}));这里先补齐movmean丢弃的两端再算残差让残差序列和原始序列对齐。movstd窗口取 14 天的原因是 Wordle 用户在周末的整体行为方差本身比工作日大7 天窗口会掩盖这种异方差14 天更稳定。阈值 2.5 不是硬规定第一轮排查先用它把最明显的干预点抓出来后面在第 4 章还要根据残差再补一轮。把第一轮抓出的日期整理成表格放进论文附录这就是模型里“干预哑变量”的来源。3. 核心模型分段趋势加周季节项的回归方程3.1 为什么选解释性强的线性回归而不是 LSTM 类模型看到“时间序列”四个字不少队伍会直接上 LSTM 或 ARIMA。C 题第一问不适合这样做理由有三点。第一它要的是“解释变化原因”论文里必须出现可以念出来的系数比如“周六比周四平均少 2.1 万人”LSTM 给不出这种结论。第二总观测数只有 300 多天把一个序列分成训练验证测试三段后深度学习模型很容易过拟合到近期异常值上。第三残差里如果仍然出现强相关的周成分说明模型结构没抓住数据生成机制换更复杂的网络也救不回来。我一般用线性模型总人数(t) 趋势(t) 周几效应(t) 干预效应(t) 噪声(t)趋势直接用分段时间基数表示段点选在 2022 年 4 月中旬附近也就是数据里总人数达到峰值的位置。这个分段点不需要特别精确后面回归系数会自适应调整。实际比赛里有人用 findchangepts 自动找分段点但 2023 C 题数据里峰值的位置肉眼可见手设段点获得的稳定性更好因为自动找点容易把 2023 年 1 月的异常低点也识别成趋势拐点。3.2 构建设计矩阵分段趋势、周哑变量与干预哑变量把上面的想法写成 MATLAB 代码先把设计矩阵的列全部构造出来n height(T); t (1:n); % 分段趋势在峰值日期 peaki 处分段 peak_date datetime(2022, 4, 15); peaki find(date peak_date, 1, first); trend1 min(t, peaki); % 上升段 trend2 max(t - peaki, 0); % 下滑段峰值之后才非零 % 周几哑变量以周一为基准生成周二到周日 6 列 dow weekday(date); % 注意 MATLAB 中 1周日 dowMat zeros(n, 6); dowMat(:,1) (dow 3); % Tue dowMat(:,2) (dow 4); % Wed dowMat(:,3) (dow 5); % Thu dowMat(:,4) (dow 6); % Fri dowMat(:,5) (dow 7); % Sat dowMat(:,6) (dow 1); % Sun % 干预哑变量outlier_idx 来自上一节 z-score 筛选 interv zeros(n, numel(outlier_idx)); for k 1:numel(outlier_idx) interv(outlier_idx(k), k) 1; end X [ones(n,1), trend1, trend2, dowMat, interv]; y double(total);这里有一个容易看晕的点MATLAB 的weekday返回 1 表示星期日和国内习惯的“周一为一周第一天”不一样。上面的代码以周一为基准周一的效应被收进截距所以后面解释的系数是“周二到周日相对于周一如何”。如果不小心把星期几映射错位周系数符号会整体乱掉。解决方式是先对dow做一次tabulate核对自己的数据映射再往下走。设计矩阵构造好后理论上可以直接用反斜杠求解最小二乘。但直接求解拿不到 t 统计量和 p 值不方便写进论文所以下一小节用fitlm封装一层。3.3 用 fitlm 求解并解读回归系数表% 构造变量名列表顺序必须和 X 列顺序一致 ivnames arrayfun((k) sprintf(Interv%d, k), ... 1:size(interv,2), UniformOutput, false); vnames [{Intercept,TrendUp,TrendDown}, ... {Tue,Wed,Thu,Fri,Sat,Sun}, ... ivnames, {Total}]; mdl fitlm(X, y, VarNames, vnames); disp(mdl.Coefficients);拟合完成后看两件事。第一件事是系数显著性pValue这一列里周系数通常全部显著因为周内波动很稳定干预系数的 p 值越大越说明这个日期被周效应吸收拉出来单独立项意义不大。第二件事是趋势系数的量级TrendUp 应当是正的、TrendDown 应当是负的如果符号反了说明 peak_date 设得太早或太晚。把主要系数的解释写进论文时可以采用类似下面这种表述变量系数含义论文里的呈现方式TrendUp峰值前上升段斜率3 月上旬至 4 月中旬日均增加约 X 人TrendDown峰值后下降段斜率之后总人数日均减少约 X 人Sat周末负效应周六比周一平均少约 X 人IntervN单个异常日增量该日总人数额外变化约 X 人提示fitlm的VarNames里最后一个名字是响应变量前面的都是预测变量。名字写错虽然不影响数值结果但会影响disp输出和后面predict的自变量名对齐所以尽量一次写对。4. 诊断与二次修正残差里漏掉的周期与异常4.1 用 autocorr 和残差图检查是否漏掉周周期第一个需要检查的是残差序列是否还存在周内规律。如果回归把周季节项吸收干净残差应该围绕 0 随机波动。一旦存在规律会表现成明显的以 7 为滞后的相关峰值。r mdl.Residuals.Raw; figure; autocorr(r, NumLags, 30);autocorr来自 Econometrics Toolbox没有这个工具箱时可以用xcorr代替但输出不如autocorr直观。重点看 7、14、21 三个滞后点如果这些点都超出 95% 置信区间说明周周期没有被完全吸收。2023 C 题的典型结果是第一次建模后 7 日滞后自相关在 0.15 左右原因是简单的“周一基准”哑变量难以处理节假日对周几年效应的平移。另一种常见情况是残差中还存在明显的正负交替这说明趋势分段点位置不准确。此时可以小幅调整 peaki 后再拟合看 AIC 是否下降。我一般不把 peaki 纳入自动参数搜索那样会让模型对噪声敏感比赛论文里也很难解释为什么参数恰好是这个值。4.2 用残差高峰补第二轮干预哑变量从残差中找出最大的几个波动点重新判断是否需要加哑变量[r_sorted, idx_sorted] sort(abs(r), descend); extra_i idx_sorted(1:3); % 每次最多加 3 个 sigma_r std(r, omitnan); new_interv zeros(n, numel(extra_i)); % 只增加值如果该日期残差大于 3 倍标准差 for k 1:numel(extra_i) if abs(r(extra_i(k))) 3 * sigma_r new_interv(extra_i(k), k) 1; end end这一步很容易陷入过度修正把残差每次最大的点都变成哑变量R² 能一直涨到 0.99但那些点没有结构也解释不了。比赛评审对这种做法容忍度很低除非能给出当时发生了什么例如某天 Wordle 单词特别难导致平均尝试次数骤升。我的原则是第一轮 z 值抓出的干预点全部保留第二轮只有在残差绝对值大于 3 倍标准差且该日期与周末或节假日相邻时才加。如果不满足宁可不加。4.3 设计矩阵不满秩时的快速排查不满秩的症状在fitlm里会显示为某些系数是 NaN或者警告Rank deficient。最常见的原因是干预哑变量所在的日期恰好是周末而周末效应已经由周哑变量解释了一部分两列高度共线。检查方法fprintf(设计矩阵秩: %d / %d\n, rank(X), size(X, 2));秩不足时不要直接删除周哑变量而是考虑把该干预日期扩展成三天窗口或者干脆移除这个干预项。2023 C 题里2022 年 11 月的黑色星期五正好落在周五从统计角度看会被周五系数吸收不需要额外设干预。多数队伍的过拟合就发生在这类操作上。5. 收尾技巧样本外验证与外推人数分布5.1 做一次最后 14 天的样本外验证第一问在论文里通常只用于解释但美赛评委常会私下补一个测试把最后 14 天数据留出来重新回归再外推。做法是把数据切成前 n-14 和最后 14 天用前段做设计矩阵拟合并预测后段nfit n - 14; mdl_fit fitlm(X(1:nfit,:), y(1:nfit), VarNames, vnames); yhat_last predict(mdl_fit, X(nfit1:end,:)); rmse_last sqrt(mean((y(nfit1:end) - yhat_last).^2)); fprintf(最后14天样本外RMSE%.1f日均总人数约%.1f\n, ... rmse_last, mean(y(nfit1:end)));如果这个值明显大于训练 RMSE 的 2 倍就说明模型在趋势拐点处不够平滑。此时考虑给 TrendDown 加一个二次项不要贸然加三次项因为三次项在外推时会把 2023 年 2 月预测成负数。作为解释类模型这个 RMSE 在 3000 到 5000 的量级已经足够。5.2 用预测结果画一张带置信带的图最终输出图是论文成败的关键。推荐用predict的置信区间参数[yhat, ci] predict(mdl, X); figure(Color,w,Position,[100 100 1000 480]); plot(date, y, k., MarkerSize, 6); hold on; plot(date, yhat, LineWidth, 1.8, Color, [0.85 0.33 0.10]); plot(date, ci(:,1), --, Color, [0.6 0.6 0.6]); plot(date, ci(:,2), --, Color, [0.6 0.6 0.6]); xlabel(日期); ylabel(总人数); legend({实际值,拟合值,置信下界,置信上界});置信带的宽度在异常日期处会明显变宽这正是干预项的体现那几天不仅偏差大方差也大。在论文里配一句“灰色虚线为 90% 置信区间干预日期附近的宽区间提示模型对突发事件的解释能力有限”就够了。5.3 给第二问留一个扩展位按比例拆总人数到各尝试次数第一问模型只输出总人数而第二问往往要求讨论不同尝试次数的人数分布。最简单的扩展是保持第一问的 Total 不变把每天某档比如 4 次猜中人数占总人数的比例做 30 天滑动平均作为基线ratio4 T.Four_Guesses ./ total; % 列名按实际数据调整 base_ratio4 movmean(ratio4, 30, Endpoints, discard);这个基线虽然粗糙但可复现、可解释。拿到第二问数据后在这个比例序列上再做周哑变量回归就比从头训练一个分类模型快得多。第一问的拟合效果直接决定了第二问比例模型的上限所以把前 5 章的结构做扎实后面所有扩展都会省时间。本文还有配套的精品资源点击获取 用MATLAB做可解释的时间序列回归——以Wordle数据为例
简介面向2023年美国大学生数学建模竞赛C题第一问提供一套基于Matlab的完整解决方案特别适合正在备赛的建模选手以及需要处理类似文本数据的科研人员。压缩包共五个文件包括两个xlsx表格分别存放原始词频数据与处理结果两个m脚本承担初始化与主流程求解任务一个mat数据文件保存了中间变量整体大小约67KB轻量易读。代码覆盖从数据读取、清洗到模型构建、结果验证与可视化的完整链路并针对Wordle数据集做了日期与词频维度的特征分析便于快速复现赛题图表。源码中附有必要的注释模块划分清晰学习者可对照赛题逐步调试。目前已有1008人学习下载说明该方案具有较高的参考价值。读者通过研读源码既能掌握Matlab处理真实题目的编程技巧也能迁移其中的数据分组统计思路至其他竞赛或科研场景中。1. 从原始数据到第一问的模型先明确这一问在解释什么2023 美赛 C 题给的是《纽约时报》Wordle 每天的玩家结果汇总列里是尝试 1 到 6 次的人数、废纸篓没猜出来的人数和总人数。第一问表面是“解释每日报告结果的变化”翻译成建模语言就是用一个数据驱动的方程把总人数随时间变化的趋势、按星期重复出现的波动以及某几天突然冒出来的异常值都解释清楚。核心不是把数字拟合得越准越好而是让模型能拆出可解释的成分。常见做法是先把数据读进来做探索性分析再用带周季节项的回归模型做主输出。MATLAB 很合适readtable 一行读 CSVdatetime 处理日期fitlm 求解很快。下面直接按能跑通的顺序写一遍源码运行环境是 MATLAB R2023b 及以上核心回归只依赖 Statistics Toolbox可选地使用 Econometrics Toolbox 画残差自相关Deep Learning Toolbox 用不到。适合准备美赛或国赛 C 题的参赛队以及想把时间序列回归做成可解释模型的工程师。2. 预处理与探索读入数据先看曲线长什么样2.1 用 readtable 和 datetime 把 Wordle 原始 CSV 读成可分析表拿到 2023 年 C 题的数据文件夹后最常见的问题是 Excel 打开看到的是乱码日期因为 CSV 里的日期是 M/d/yyyy 格式。用 MATLAB 读的时候不要用 importdata建议直接 readtable 然后自己控制日期格式clear; rng(42); % 读取原始数据文件路径改成你本机实际路径 T readtable(Wordle_Data_March_2022_to_Jan_2023.csv, ... TextType, string); % 日期列为 M/d/yyyy 格式要显式指定否则整列是 NaT T.Date datetime(T.Date, InputFormat, M/d/yyyy); T sortrows(T, Date); % 取当天总人数列原数据里这一列名通常叫 Total total T.Total_Number_of_Reported_Results; date T.Date; fprintf(样本范围: %s 到 %s, 共 %d 天\n, ... datestr(min(date)), datestr(max(date)), height(T));代码里先用TextType让字符串列以 string 类型进入内存避免后面字符数组拼接出问题datetime的InputFormat要和 CSV 里实际格式严格一致不然会得到一整列 NaT。很多参赛队在这里翻车Excel 把日期显示成 3/1/2022MATLAB 默认格式是 yyyy-MM-dd不加 InputFormat 会全部解析失败。读进来之后建议顺手做两层检查sum(isnat(date))必须为 0sum(isnan(total))必须为 0。这两个数只要有一个非零后面回归结果就不该继续往下走得先回去看是原始文件空行还是合并单元格引起的。2.2 用 7 天滑动平均将趋势和周内波动分离开原始总人数序列是一个高频时间序列直接看曲线只能看到“2022 年 4 月高、之后缓慢下降、偶尔有尖峰”看不出周末效应。我一般会先画一张原始曲线和 7 天滑动平均叠在一起的图% movmean 用 7 天窗口等价于消除周内周期 trend7 movmean(total, 7, Endpoints, discard); figure(Color, w, Position, [100 100 900 420]); plot(date, total, Color, [0.75 0.75 0.75]); hold on; plot(date(4:end-3), trend7, LineWidth, 2, ... Color, [0.85 0.33 0.10]); xlabel(日期); ylabel(当天报告总人数); legend({原始总人数, 7天滑动平均}, Location, best); grid on; box on;movmean的Endpoints指定两端点如何处理discard会让输出首尾各少 3 个观测。这样得到的 trend7 去掉的是以 7 天为周期的波动剩下的变化就主要是长期趋势加突发事件。对 2023 年 C 题的数据你会观察到日期区间原始现象建模时对应处理2022 年 3 月到 4 月中旬总人数快速上升分段上升趋势2022 年 4 月到 2023 年 1 月总人数稳步下滑分段下滑趋势2023 年 1 月初几天明显低点干预哑变量某些星期一或星期二出现低点周系数取负值这个表可以作为预处理的第一个输出放进论文评委看到的不只是一张图还能看出你已经把“趋势 周效应 异常”三个成分分开。2.3 对偏差过大的日期做 z-score 检查并记录干预点滑动平均能看趋势但没法客观地标出哪些日期算异常。更稳妥的做法是先算残差原始值减 trend7再用滚动标准差做 z-score。具体实现% 先用 NaN 把 movmean 丢弃的两端补回来再做逐点残差 resid total - [NaN(3,1); trend7; NaN(3,1)]; sd14 movstd(total, 14, Endpoints, discard); z abs(resid(4:end-3) ./ sd14); % 找出 z 超过 2.5 的日期 outlier_idx find(z 2.5); outlier_dates date(4:end-3); disp(table(outlier_dates(outlier_idx), total(4:end-3), z(outlier_idx), ... VariableNames, {日期, 总人数, z值}));这里先补齐movmean丢弃的两端再算残差让残差序列和原始序列对齐。movstd窗口取 14 天的原因是 Wordle 用户在周末的整体行为方差本身比工作日大7 天窗口会掩盖这种异方差14 天更稳定。阈值 2.5 不是硬规定第一轮排查先用它把最明显的干预点抓出来后面在第 4 章还要根据残差再补一轮。把第一轮抓出的日期整理成表格放进论文附录这就是模型里“干预哑变量”的来源。3. 核心模型分段趋势加周季节项的回归方程3.1 为什么选解释性强的线性回归而不是 LSTM 类模型看到“时间序列”四个字不少队伍会直接上 LSTM 或 ARIMA。C 题第一问不适合这样做理由有三点。第一它要的是“解释变化原因”论文里必须出现可以念出来的系数比如“周六比周四平均少 2.1 万人”LSTM 给不出这种结论。第二总观测数只有 300 多天把一个序列分成训练验证测试三段后深度学习模型很容易过拟合到近期异常值上。第三残差里如果仍然出现强相关的周成分说明模型结构没抓住数据生成机制换更复杂的网络也救不回来。我一般用线性模型总人数(t) 趋势(t) 周几效应(t) 干预效应(t) 噪声(t)趋势直接用分段时间基数表示段点选在 2022 年 4 月中旬附近也就是数据里总人数达到峰值的位置。这个分段点不需要特别精确后面回归系数会自适应调整。实际比赛里有人用 findchangepts 自动找分段点但 2023 C 题数据里峰值的位置肉眼可见手设段点获得的稳定性更好因为自动找点容易把 2023 年 1 月的异常低点也识别成趋势拐点。3.2 构建设计矩阵分段趋势、周哑变量与干预哑变量把上面的想法写成 MATLAB 代码先把设计矩阵的列全部构造出来n height(T); t (1:n); % 分段趋势在峰值日期 peaki 处分段 peak_date datetime(2022, 4, 15); peaki find(date peak_date, 1, first); trend1 min(t, peaki); % 上升段 trend2 max(t - peaki, 0); % 下滑段峰值之后才非零 % 周几哑变量以周一为基准生成周二到周日 6 列 dow weekday(date); % 注意 MATLAB 中 1周日 dowMat zeros(n, 6); dowMat(:,1) (dow 3); % Tue dowMat(:,2) (dow 4); % Wed dowMat(:,3) (dow 5); % Thu dowMat(:,4) (dow 6); % Fri dowMat(:,5) (dow 7); % Sat dowMat(:,6) (dow 1); % Sun % 干预哑变量outlier_idx 来自上一节 z-score 筛选 interv zeros(n, numel(outlier_idx)); for k 1:numel(outlier_idx) interv(outlier_idx(k), k) 1; end X [ones(n,1), trend1, trend2, dowMat, interv]; y double(total);这里有一个容易看晕的点MATLAB 的weekday返回 1 表示星期日和国内习惯的“周一为一周第一天”不一样。上面的代码以周一为基准周一的效应被收进截距所以后面解释的系数是“周二到周日相对于周一如何”。如果不小心把星期几映射错位周系数符号会整体乱掉。解决方式是先对dow做一次tabulate核对自己的数据映射再往下走。设计矩阵构造好后理论上可以直接用反斜杠求解最小二乘。但直接求解拿不到 t 统计量和 p 值不方便写进论文所以下一小节用fitlm封装一层。3.3 用 fitlm 求解并解读回归系数表% 构造变量名列表顺序必须和 X 列顺序一致 ivnames arrayfun((k) sprintf(Interv%d, k), ... 1:size(interv,2), UniformOutput, false); vnames [{Intercept,TrendUp,TrendDown}, ... {Tue,Wed,Thu,Fri,Sat,Sun}, ... ivnames, {Total}]; mdl fitlm(X, y, VarNames, vnames); disp(mdl.Coefficients);拟合完成后看两件事。第一件事是系数显著性pValue这一列里周系数通常全部显著因为周内波动很稳定干预系数的 p 值越大越说明这个日期被周效应吸收拉出来单独立项意义不大。第二件事是趋势系数的量级TrendUp 应当是正的、TrendDown 应当是负的如果符号反了说明 peak_date 设得太早或太晚。把主要系数的解释写进论文时可以采用类似下面这种表述变量系数含义论文里的呈现方式TrendUp峰值前上升段斜率3 月上旬至 4 月中旬日均增加约 X 人TrendDown峰值后下降段斜率之后总人数日均减少约 X 人Sat周末负效应周六比周一平均少约 X 人IntervN单个异常日增量该日总人数额外变化约 X 人提示fitlm的VarNames里最后一个名字是响应变量前面的都是预测变量。名字写错虽然不影响数值结果但会影响disp输出和后面predict的自变量名对齐所以尽量一次写对。4. 诊断与二次修正残差里漏掉的周期与异常4.1 用 autocorr 和残差图检查是否漏掉周周期第一个需要检查的是残差序列是否还存在周内规律。如果回归把周季节项吸收干净残差应该围绕 0 随机波动。一旦存在规律会表现成明显的以 7 为滞后的相关峰值。r mdl.Residuals.Raw; figure; autocorr(r, NumLags, 30);autocorr来自 Econometrics Toolbox没有这个工具箱时可以用xcorr代替但输出不如autocorr直观。重点看 7、14、21 三个滞后点如果这些点都超出 95% 置信区间说明周周期没有被完全吸收。2023 C 题的典型结果是第一次建模后 7 日滞后自相关在 0.15 左右原因是简单的“周一基准”哑变量难以处理节假日对周几年效应的平移。另一种常见情况是残差中还存在明显的正负交替这说明趋势分段点位置不准确。此时可以小幅调整 peaki 后再拟合看 AIC 是否下降。我一般不把 peaki 纳入自动参数搜索那样会让模型对噪声敏感比赛论文里也很难解释为什么参数恰好是这个值。4.2 用残差高峰补第二轮干预哑变量从残差中找出最大的几个波动点重新判断是否需要加哑变量[r_sorted, idx_sorted] sort(abs(r), descend); extra_i idx_sorted(1:3); % 每次最多加 3 个 sigma_r std(r, omitnan); new_interv zeros(n, numel(extra_i)); % 只增加值如果该日期残差大于 3 倍标准差 for k 1:numel(extra_i) if abs(r(extra_i(k))) 3 * sigma_r new_interv(extra_i(k), k) 1; end end这一步很容易陷入过度修正把残差每次最大的点都变成哑变量R² 能一直涨到 0.99但那些点没有结构也解释不了。比赛评审对这种做法容忍度很低除非能给出当时发生了什么例如某天 Wordle 单词特别难导致平均尝试次数骤升。我的原则是第一轮 z 值抓出的干预点全部保留第二轮只有在残差绝对值大于 3 倍标准差且该日期与周末或节假日相邻时才加。如果不满足宁可不加。4.3 设计矩阵不满秩时的快速排查不满秩的症状在fitlm里会显示为某些系数是 NaN或者警告Rank deficient。最常见的原因是干预哑变量所在的日期恰好是周末而周末效应已经由周哑变量解释了一部分两列高度共线。检查方法fprintf(设计矩阵秩: %d / %d\n, rank(X), size(X, 2));秩不足时不要直接删除周哑变量而是考虑把该干预日期扩展成三天窗口或者干脆移除这个干预项。2023 C 题里2022 年 11 月的黑色星期五正好落在周五从统计角度看会被周五系数吸收不需要额外设干预。多数队伍的过拟合就发生在这类操作上。5. 收尾技巧样本外验证与外推人数分布5.1 做一次最后 14 天的样本外验证第一问在论文里通常只用于解释但美赛评委常会私下补一个测试把最后 14 天数据留出来重新回归再外推。做法是把数据切成前 n-14 和最后 14 天用前段做设计矩阵拟合并预测后段nfit n - 14; mdl_fit fitlm(X(1:nfit,:), y(1:nfit), VarNames, vnames); yhat_last predict(mdl_fit, X(nfit1:end,:)); rmse_last sqrt(mean((y(nfit1:end) - yhat_last).^2)); fprintf(最后14天样本外RMSE%.1f日均总人数约%.1f\n, ... rmse_last, mean(y(nfit1:end)));如果这个值明显大于训练 RMSE 的 2 倍就说明模型在趋势拐点处不够平滑。此时考虑给 TrendDown 加一个二次项不要贸然加三次项因为三次项在外推时会把 2023 年 2 月预测成负数。作为解释类模型这个 RMSE 在 3000 到 5000 的量级已经足够。5.2 用预测结果画一张带置信带的图最终输出图是论文成败的关键。推荐用predict的置信区间参数[yhat, ci] predict(mdl, X); figure(Color,w,Position,[100 100 1000 480]); plot(date, y, k., MarkerSize, 6); hold on; plot(date, yhat, LineWidth, 1.8, Color, [0.85 0.33 0.10]); plot(date, ci(:,1), --, Color, [0.6 0.6 0.6]); plot(date, ci(:,2), --, Color, [0.6 0.6 0.6]); xlabel(日期); ylabel(总人数); legend({实际值,拟合值,置信下界,置信上界});置信带的宽度在异常日期处会明显变宽这正是干预项的体现那几天不仅偏差大方差也大。在论文里配一句“灰色虚线为 90% 置信区间干预日期附近的宽区间提示模型对突发事件的解释能力有限”就够了。5.3 给第二问留一个扩展位按比例拆总人数到各尝试次数第一问模型只输出总人数而第二问往往要求讨论不同尝试次数的人数分布。最简单的扩展是保持第一问的 Total 不变把每天某档比如 4 次猜中人数占总人数的比例做 30 天滑动平均作为基线ratio4 T.Four_Guesses ./ total; % 列名按实际数据调整 base_ratio4 movmean(ratio4, 30, Endpoints, discard);这个基线虽然粗糙但可复现、可解释。拿到第二问数据后在这个比例序列上再做周哑变量回归就比从头训练一个分类模型快得多。第一问的拟合效果直接决定了第二问比例模型的上限所以把前 5 章的结构做扎实后面所有扩展都会省时间。本文还有配套的精品资源点击获取