简介这份2023年美国大学生数学建模竞赛C题的M奖Meritorious论文围绕Wordle游戏数据展开完整的建模与预测分析适合备战美赛的选手、需要参考英文论文写作规范的学生以及做数据分析课程设计或毕业设计的读者使用。全文为英文覆盖数据清洗与缺失值处理、9项单词属性的特征提取与Spearman相关性分析、时间序列区间预测、偏最小二乘回归PLS建模、K-means聚类结合肘部法则与轮廓系数确定词性与难度分组等环节并给出MAE在34、MSE在2224之间的误差评估与混淆矩阵热图还以单词EERIE为例算出4.635的难度得分并归入第五难度区间展示从聚类到评分的完整实现路径文末附有模型不足与不确定性讨论可作为完整赛题方案、论文结构与英文表达的学习范本。压缩包内共1个PDF文件约2.18MB为论文全文便于直接阅读与打印。目前已有170人学习下载适合希望快速理解M奖论文行文逻辑与建模思路的参赛者参考。1. 一份 25 页的 C 题论文卡住复现的往往不是公式去年有人拿着一份建模赛 C 题 M 奖论文来找我说只想把「3 月 1 日的结果数区间」这一问跑出来公式看得懂代码却怎么也复现不出[17363, 18994]这个区间。翻到数据那一节才发现卡点根本不在 ARIMA而在原始表里两个只有四个字母的单词和两处断崖式异常值上。这份论文做的是某每日五字母猜词游戏的公开数据用历史数据做时间序列区间预测、用九个词属性解释困难模式占比、用偏最小二乘预测某个词未来七档尝试次数的分布、再用 K-means 把单词按难度分级最后还留了一节讲数据里的有趣特征。它特别适合两类人一类是拿它当毕业设计或课程设计骨架、想找一条「清洗—特征—时序—回归—聚类—评估」完整链路的学生另一类是想看一份小样本多因变量建模怎么落地、又不想从零搭工程的开发者。下面按我拆它的顺序讲先看骨架再逐模块跑。2. 论文骨架拆解25 页目录怎么对应到一条可执行的数据流拿到这种论文我第一件事不是读模型是把目录抄成一张「章节 → 输入 → 输出 → 落地产物」的表。论文的章节编号本身已经暴露了作者的设计意图数据预处理单独成章放在最前面说明后四个问题共用同一份清洗结果第四问有趣特征放在模型之后说明它是对前面建模过程的观察沉淀不是独立任务。把这张表列出来你就知道哪一段该写脚本、哪一段只是结论陈述。论文章节对应任务输入输出II. Data preprocessing清洗 可视化359 行原始表12 列干净表、趋势图III. Problem 1时序预测 属性相关时间序列、9 个词属性3/1 预测区间、Spearman 系数IV. Problem 2多因变量回归词属性矩阵、七档百分比每档百分比的预测值、MAE/MSEV. Problem 3无监督难度分级词属性 预测百分比5 档难度、EERIE 得分、混淆矩阵VI. Problem 4特征观察全量数据若干规律性结论从数据流角度看整篇论文其实就一条主干一张宽表word, date, num_reported, num_hard, p1..p6, pX先被清洗成稳定序列再被拆成两个方向——时间维度走 ARIMA单词维度走属性提取 PLS最后两条线在第三问汇合用 PLS 预测出的七档百分比当权重算难度分再用 K-means 分档验证。项目如果要做成可复跑的代码仓我一般会按这个主干拆成四个文件wordle-c2023/ ├── data/ │ ├── raw/ # 原始 359 行数据只读不改 │ └── clean/ # 清洗后的中间结果可随时重建 ├── src/ │ ├── preprocess.py # 第二章纠错、异常值、百分比归一化 │ ├── forecast.py # 第三章ARIMA 定阶与区间预测 │ ├── features.py # 第三章9 个词属性 Spearman │ └── difficulty.py # 第四、五章PLS K-means 混淆矩阵 └── outputs/ ├── figures/ # 论文里的曲线图、饼图、散点图 └── tables/ # 预测区间、难度分、混淆矩阵这个结构看起来朴素但有两个好处raw目录只读保证任何时候都能从原始数据重跑一遍每一章对应一个脚本改模型不会牵一发动全身。论文里没有给出代码组织方式这是按常规工程习惯补的你完全可以把forecast.py和features.py合成一个只要保证清洗脚本先跑就行。2.1 第一刀清洗两个四字母单词和两处断崖式异常论文里写得很轻描淡写「发现有两个单词只有四个字母查证后修正为 clean 和 trash」。真跑起来这一步不能省因为后面所有基于「五字母」假设提取的属性首字母、尾字母、字母位置分布都会因为长度不对而错位。这类错误在公开数据集里很常见通常是录入时漏了字母处理方式就是人工核对后做一张映射表别用程序去猜。import pandas as pd import numpy as np df pd.read_csv(data/raw/wordle.csv) df[date] pd.to_datetime(df[date]) df df.sort_values(date).reset_index(dropTrue) # 1) 长度异常修正只对确认过的两个词做映射不做模糊匹配 typo_map {clen: clean, tash: trash} df[word] df[word].str.lower().str.strip().replace(typo_map) # 修正后必须全表满足五字母否则后面属性提取会静默出错 bad df[df[word].str.len() ! 5] assert bad.empty, f仍有非五字母单词{bad[word].tolist()} # 2) 断崖式异常用前后各一个点的均值替换 def fix_spike(frame, col, ts, window1): idx frame.index[frame[date] pd.Timestamp(ts)] if len(idx) 0: raise KeyError(f{ts} 不在数据集中) i idx[0] left, right frame.loc[i - window, col], frame.loc[i window, col] frame.loc[i, col] int(round((left right) / 2)) return frame df fix_spike(df, num_reported, 2022-11-30) # 2569 - (2373922628)/2 df fix_spike(df, num_hard, 2022-02-13) # 3249 - (931010343)/2这段代码有三个参数值得说清楚。typo_map必须是显式白名单我见过有人用「长度不等于 5 就丢弃」的写法结果把两个本来能救回来的样本扔了359 行变成 357 行后面的时序长度就和论文对不上。window1表示只取左右各一个邻居因为这两处异常都是单点尖刺窗口开大会把正常波动也抹平。round之后再取整是因为论文里的 23184 和 9827 都是四舍五入的结果不取整会在后续画图时出现小数刻度和论文图对不上。异常值为什么不能直接删因为这是时间序列删掉一行等于把 11 月 30 日和 12 月 1 日强行接在一起ARIMA 的滞后项会因此错位一天。均值替换虽然粗糙但保住了时间轴的连续性这是这套论文在数据量只有 359 行时比较务实的选择。2.2 百分比归一化98% 和 102% 都得变成 100%第二个坑更隐蔽。论文提到七档百分比之和在数据集中出现了 98%、99%、100%、101%、102% 五种情况。如果你直接拿原始百分比做 PLS 的因变量模型会学到一种假信号某些日期的「总概率质量」偏大或偏小但这个偏差其实来自四舍五入和平台统计口径不是玩家行为差异。pct_cols [p1, p2, p3, p4, p5, p6, pX] # 每行七档百分比求和再逐列除以总和强制归一 row_sum df[pct_cols].sum(axis1) df[pct_cols] df[pct_cols].div(row_sum, axis0) * 100 # 归一后校验每行之和与 100 的偏差应小于 1e-6 err (df[pct_cols].sum(axis1) - 100).abs().max() assert err 1e-6, f归一化未收敛最大偏差 {err:.6f} # 困难模式占比用人数直接相除比用百分比相减更稳 df[hard_ratio] df[num_hard] / df[num_reported] * 100这里div(row_sum, axis0)是按行广播axis0表示除数是每一行的标量容易写反成axis1写反了就是按列归一结果是每一列加起来等于 100看起来也「像对的」但完全破坏了样本间的比例关系属于典型的静默错误。断言里的1e-6是浮点容差太小会因为二进制表示误差误报太大又查不出问题。困难模式占比我没有用p1..pX那七列去推而是直接用num_hard / num_reported原因是论文第一问要分析的就是「困难模式得分百分比」和词属性的关系用人数相除得到的比例不受归一化舍入影响。这一步做完后面第三问的 Spearman 相关才有意义。2.3 可视化初检三张图其实是在给模型做体检论文第二章给了三组图七档百分比堆叠图、未成功挑战人数散点图、困难模式通关占比趋势图。很多人把这三张图当装饰其实它们各自对应一个建模假设的检验。堆叠图里黄色区域第四次尝试成功占比最大说明分布的重心在 3 到 5 次之间这直接决定了后面 PLS 的因变量不能当成七个独立回归硬拟合而要接受它们之间存在此消彼长的约束。未成功散点图里绝大多数日期低于 5%只有几十天超过 10%说明数据里存在少量「难词日」这些点在 K-means 分档时会被推到高档位如果直接当噪声剔掉难度分档的边界会整体左移。困难模式占比趋势图是单调上升的这一条最关键它意味着hard_ratio不是一个平稳序列如果你把它当因变量做时序回归必须先去趋势否则模型会把「时间」误当成「词属性」的解释力。我一般会在清洗脚本最后加一段快速出图只用来肉眼确认不追求论文里的排版import matplotlib.pyplot as plt fig, axes plt.subplots(1, 3, figsize(15, 4)) axes[0].stackplot(df[date], [df[c] for c in pct_cols]) axes[0].set_title(七档尝试次数占比归一后) axes[1].scatter(df[date], df[pX], s8) axes[1].axhline(5, ls--, cr) # 5% 参考线 axes[1].set_title(未成功挑战占比) axes[2].plot(df[date], df[hard_ratio]) axes[2].set_title(困难模式占比趋势) plt.tight_layout() plt.savefig(outputs/figures/preprocess_check.png, dpi150)axhline(5)那条红线不是论文原图里的是我自己加的参考线用来快速判断有多少天越过了 5% 这个量级。三张图放在一张画布上跑完清洗脚本顺手看一眼比等到建模阶段报错再回头查要省时间得多。3. 建模主干ARIMA 区间预测与 PLS 多因变量回归怎么落地清洗完之后论文的建模部分其实分成两条互不干扰的线。一条是时间维度只关心num_reported这个单变量序列目标是给出 2023 年 3 月 1 日的区间另一条是单词维度把每个词压成九个属性去解释困难模式占比再去预测七档百分比。两条线在代码上应该完全解耦因为它们的样本单位不一样前者是「天」后者是「词」混在一个脚本里最容易出现的错误就是把日期索引和单词索引对错行。3.1 ARIMA 定阶为什么最后落在低阶 (p,d,q)论文用 ARIMA(p,d,q) 做区间预测给出的结果是 [17363, 18994]。这三个参数里d 是差分阶数作用是让序列平稳p 是自回归项数q 是滑动平均项数。论文正文推导了(1-ΣφᵢLⁱ)(1-L)^d Xₜ (1ΣθᵢLⁱ)εₜ这个表达式落到代码里其实就是一次网格搜索加 AIC 选优。from statsmodels.tsa.arima.model import ARIMA # 以日期为索引缺失日期补全后再插值避免时间轴断裂 y df.set_index(date)[num_reported].asfreq(D).interpolate() best_aic, best_order, best_model np.inf, None, None for p in range(0, 4): for d in range(0, 3): for q in range(0, 4): try: m ARIMA(y, order(p, d, q), enforce_stationarityFalse, enforce_invertibilityFalse).fit() except Exception: continue # 不收敛的组合直接跳过 if m.aic best_aic: best_aic, best_order, best_model m.aic, (p, d, q), m print(最优阶数, best_order, AIC, round(best_aic, 2))参数上enforce_stationarityFalse和enforce_invertibilityFalse是为了让搜索空间更宽代价是可能选出物理意义不强的参数组合所以实战里我会把搜索上限压在 p4、q4。asfreq(D)加interpolate()这两步不能省因为 359 天里存在平台未记录的空档索引不连续时 statsmodels 的滞后计算会按行号而不是按日期走预测出来的区间会偏窄。拿到最优模型后预测区间这样出forecast best_model.get_forecast(steps(pd.Timestamp(2023-03-01) - y.index[-1]).days) pred forecast.predicted_mean.loc[2023-03-01] ci forecast.conf_int(alpha0.05).loc[2023-03-01] # 95% 置信区间 print(f点预测 {pred:.0f}区间 [{ci.iloc[0]:.0f}, {ci.iloc[1]:.0f}])alpha0.05对应 95% 置信水平论文里的 [17363, 18994] 就是这个口径下的结果。这里有个容易混淆的点置信区间和预测区间不是一回事前者描述的是「均值的估计误差」后者还要叠加单点观测的随机波动宽度更大。论文在正文里两个都建了代码上分别对应get_forecast().conf_int()和get_prediction().conf_int()如果你复现出来的区间比论文窄一截先检查是不是只取了均值置信区间。3.2 九个词属性怎么挑Spearman 比 Pearson 更适合这个场景论文提取了 9 个独立变量对困难模式得分百分比做 Spearman 相关分析结论是「这些属性有一定影响」。九个属性论文没有全部列出名称但从上下文能推断出大致方向词长这里恒为 5可作占位、重复字母个数、元音个数、字母频率总和、首字母频率、尾字母频率、是否含常见前缀、是否含常见后缀、字母分布熵。这套属性在猜词类游戏分析里是常规做法因为玩家猜词的难度主要来自「字母是否常见」和「结构是否有提示」。为什么用 Spearman 而不是 Pearson因为七档百分比是有界的且分布明显右偏第四次尝试那一档占比最大Pearson 假设线性关系和近似正态在这种有界比例数据上容易被极端值带偏。Spearman 只看秩对单调非线性更稳。from scipy.stats import spearmanr feat_cols [dup_cnt, vowel_cnt, freq_sum, first_freq, last_freq, has_prefix, has_suffix, entropy, uniq_cnt] rows [] for col in feat_cols: rho, pval spearmanr(df[col], df[hard_ratio]) rows.append({attr: col, rho: round(rho, 3), p: round(pval, 4)}) corr pd.DataFrame(rows).sort_values(rho, keyabs, ascendingFalse) print(corr)rho的绝对值超过 0.3 基本可以认为有中等程度关联p小于 0.05 说明这个关联不太可能是偶然。这里要注意九个属性彼此之间可能高度相关比如vowel_cnt和entropy单变量相关分析只能筛出候选不能直接下「因果」结论。论文的表述是「some effect」这个措辞是准确的因为后面 PLS 用的是全部九个属性一起进模型靠潜变量去吸收共线性而不是靠单变量筛选。我实际跑这类数据时会在相关分析后加一步 VIF 检查如果某个属性 VIF 超过 10就在 PLS 里降低它的权重或者干脆剔除否则潜变量个数会虚高。3.3 PLS 回归样本 359 行、七个因变量普通最小二乘为什么不够用第二问要把七个百分比1 到 6 次成功外加 X 未成功同时作为因变量预测自变量是九个词属性。这是一个典型的多因变量小样本问题样本量 359自变量 9 个因变量 7 个而且七个因变量之和恒等于 100存在完全共线性——你不可能同时独立地预测七个数因为它们被一个约束绑死了。普通多元线性回归在这里会直接翻车系数矩阵不满秩。论文选 PLS偏最小二乘是对的。PLS 的思路是先对 X 和 Y 同时做潜变量提取类似 PCA但提取方向要最大化与 Y 的协方差再在潜变量空间里做回归。它把主成分分析、多元线性回归和误差分析揉在一起对样本量适中、变量间共线性强的数据特别友好。from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import cross_val_predict, KFold from sklearn.metrics import mean_absolute_error, mean_squared_error X df[feat_cols].values Y df[pct_cols].values # 5 折交叉验证选潜变量个数避免用测试集调参 kf KFold(n_splits5, shuffleTrue, random_state42) best_n, best_mae 1, np.inf for n in range(1, 8): pls PLSRegression(n_componentsn, scaleTrue) Y_pred cross_val_predict(pls, X, Y, cvkf) mae mean_absolute_error(Y, Y_pred) if mae best_mae: best_n, best_mae n, mae final PLSRegression(n_componentsbest_n, scaleTrue).fit(X, Y) Y_hat final.predict(X) print(潜变量个数, best_n, MAE, round(mean_absolute_error(Y, Y_hat), 2), MSE, round(mean_squared_error(Y, Y_hat), 2))scaleTrue是关键参数它会对 X 和 Y 都做标准化。九个属性的量纲差别很大——freq_sum可能是几十has_prefix只有 0/1——不标准化的话量纲大的属性会独占第一个潜变量。n_components论文没有明说取了几个但常规做法是在 1 到 5 之间用交叉验证选选太大就会过拟合交叉验证 MAE 会先降后升。论文给出的 MAE 在 3 到 4 之间、MSE 在 22 到 24 之间这个量级是合理的七档百分比的均值大概在 14 左右MAE 3 到 4 意味着平均每档误差三到四个百分点。要注意 MAE 和 MSE 的单位是「百分点」不是比例所以 22 到 24 的 MSE 开根号后 RMSE 大约 4.7和 MAE 同量级说明误差分布没有特别重的长尾。4. 难度分级与验证K-means 分档、加权打分和 0.941 的准确率第三问是整个论文里最容易写得好看、也最容易复现翻车的一问。任务本身不复杂单词没有现成的难度标签所以用无监督学习分档分完档之后用 PLS 预测的七档百分比算一个加权难度分最后用混淆矩阵验证分档的准确率。三个环节环环相扣任何一环的标签顺序错了最后那个 0.941 就出不来。4.1 聚类数怎么定肘部法则看拐点轮廓系数看分离度K-means 最大的问题是 k 要自己给。论文用了两个指标交叉验证肘部法则和轮廓系数。肘部法则看的是簇内平方和inertia随 k 变化的曲线拐点处就是「再加一个簇带来的收益开始变小」的位置轮廓系数衡量的是每个点到本簇的紧密度和到最近邻簇的分离度取值 -1 到 1越高越好。from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler from sklearn.metrics import silhouette_score Xw StandardScaler().fit_transform(df[feat_cols]) # 聚类前必须标准化 inertia, sil [], {} for k in range(2, 11): km KMeans(n_clustersk, n_init20, random_state42).fit(Xw) inertia.append(km.inertia_) sil[k] silhouette_score(Xw, km.labels_) # 肘部相邻两点斜率下降最剧烈的位置 drops np.diff(inertia) / inertia[:-1] elbow_k int(np.argmin(drops)) 2 # 2 是因为 k 从 2 开始 print(肘部候选 k , elbow_k, 轮廓系数 , round(sil[elbow_k], 3))n_init20表示用 20 组不同初始质心各跑一次取最优K-means 对初始化敏感默认的 10 次在小样本上偶尔会掉进局部最优。random_state42是为了可复现换了随机种子簇编号会变但分档结构一般不变——这一点很重要因为后面混淆矩阵的标签对齐依赖簇编号固定种子能省掉一堆麻烦。论文最后定的是五档难度这个数字通常不是肘部法则直接给出的而是肘部候选常见是 3 到 4和轮廓系数峰值之间折中再结合「五档在业务上好解释」这个约束定下来的。我一般会把 inertia 和轮廓系数的曲线都画出来放在一起看两个指标都支持的那个 k 才靠谱。4.2 难度分怎么算EERIE 的 4.635 是加权期望分完档只是有了「组」组内怎么排序、怎么给出一个连续分数论文用的是加权公式把 PLS 预测出的七档百分比当作权重乘上对应的尝试次数求和得到一个期望尝试次数值越大表示越难。# p1..p6 对应 1~6 次成功pX 是未成功权重记为 7 weights np.array([1, 2, 3, 4, 5, 6, 7]) def difficulty_score(pct_row): pct np.asarray(pct_row, dtypefloat) pct pct / pct.sum() # 防止归一化残差 return float((pct * weights).sum()) eer df.loc[df[word] eerie, pct_cols].values[0] print(EERIE 难度分 , round(difficulty_score(eer), 3)) # 论文中为 4.635权重 7 是给未成功那一档的这个取值不是唯一的取 7 相当于假设「失败比第六次成功再难一档」你也可以取 6.5 或 8难度分的绝对值会变但排序基本稳定。EERIE 这个词有重复字母、元音密集属于典型的「看着简单猜着难」4.635 落在五档里的高档位和论文结论一致。这一步有个容易忽略的细节PLS 预测出来的百分比可能为负或者和不为 100所以pct / pct.sum()这行归一化不能省否则加权和会出现负数或者超过 7 的异常值。4.3 混淆矩阵 0.941验证的不是聚类是「预测能不能落回原档」论文用混淆矩阵热图评估模型准确率 0.941。这里要理解清楚它在验证什么K-means 本身是无监督的没有真标签所以混淆矩阵的标签来源只能是「用预测百分比重新算难度分再按难度区间反推档位」和原始聚类结果对比。换句话说它验证的是「PLS 预测 → 加权打分 → 分档」这条链路和直接聚类的结果是否一致一致性高说明 PLS 的预测没有把词的结构信息丢掉。from sklearn.metrics import confusion_matrix, accuracy_score # 用全量数据重新算难度分并按分位切五档 scores np.array([difficulty_score(row) for row in df[pct_cols].values]) bins np.quantile(scores, [0, 0.2, 0.4, 0.6, 0.8, 1.0]) pred_level np.digitize(scores, bins[1:-1]) 1 # true_level 来自 K-means 的簇编号按簇内均分排序后映射到 1~5 cm confusion_matrix(true_level, pred_level, labels[1, 2, 3, 4, 5]) acc accuracy_score(true_level, pred_level) print(cm) print(准确率, round(acc, 3)) # 论文中为 0.941np.quantile用等频分箱保证五档样本量接近这比等宽分箱更稳因为难度分是右偏的等宽分箱会把大部分词塞进低档。labels[1,2,3,4,5]显式指定标签顺序防止混淆矩阵行列顺序错位——这是复现这类结果最常见的翻车点代码跑通但数值对不上八成是标签顺序反了。准确率 0.941 意味着 359 个词里大约 21 个被分到了相邻档位这个误差水平在无监督分档里算相当好。看混淆矩阵时重点看对角线附近如果误差主要集中在 2 和 3 之间说明这两档的边界本来就模糊属于正常如果误差散落在对角线上那就要怀疑属性提取有问题。4.4 第四问的特征观察几个值得单独记下来的规律第四问没有建模是对全量数据的观察总结论文列了四条时间上的相关性、首字母频率、尝试次数的近似正态分布、以及某些特殊词缀的异常。这几条其实都能从前面清洗后的数据里直接验证。时间相关性指的是结果数序列存在周内波动周末参与人数偏高这也是 ARIMA 里 d 取小的原因之一过度差分会把周效应也差掉。首字母频率说的是以常见字母开头的词玩家第一次猜中的概率更高这会直接影响p1这一档的百分比分布。尝试次数近似正态这一点和第二章的堆叠图对得上重心在第四次。特殊词缀异常是指含常见后缀的词反而更难因为玩家容易先入为主地往错误方向猜。这四条结论的价值不在于多深奥而在于它们都是「建模过程中顺手能看到」的副产物。我拆这类论文时习惯把第四问当成检查清单如果你复现出来的数据也能观察到这四条说明前面的清洗和特征提取方向是对的如果一条都看不到先回头查归一化和异常值处理。5. 避坑与排查复现这套代码时最容易翻车的五处这一章是我自己踩过的坑按「现象 → 原因 → 解决」写。这五个问题有一个共同特点代码不报错结果看起来正常但数值就是和论文对不上。5.1 数据层的三个坑坑一百分比归一化写成了按列归一。现象是模型能跑MAE 也不大但七档百分比的预测值加起来不等于 100强行解释成「模型误差」。原因是df[pct_cols].div(row_sum, axis0)里的axis写反了写成了axis1变成每个列向量除以一个标量序列实际效果是按列缩放。解决方式很简单在归一化后加一行断言(df[pct_cols].sum(axis1) - 100).abs().max() 1e-6让它在清洗阶段就炸出来而不是等到建模阶段猜。坑二异常值替换时用错了邻居。现象是替换后序列在异常点附近出现一个明显的台阶ARIMA 的预测区间整体偏移。原因是代码里用了df.loc[i-1]和df.loc[i1]取邻居但这个 DataFrame 的索引在排序后没有重置loc取的是标签不是位置取到了别的日期。解决方式是清洗开始时统一reset_index(dropTrue)替换函数里改用iloc[i-window]和iloc[iwindow]并且把替换前后的值打印出来核对2569 变 23184、3249 变 9827这两个数对上了再往下走。坑三时间序列没补全就进 ARIMA。现象是模型能拟合但预测到 3 月 1 日的步数和论文差几天区间宽度明显偏窄。原因是原始数据里有些日期平台没记录直接set_index(date)后索引不连续statsmodels 按行号推算滞后把不连续的日子当成了连续。解决方式是asfreq(D)补齐日历日再interpolate()插值插值方法用线性就够因为缺的是零星单日不是整段缺失。5.2 模型层的两个坑坑四PLS 忘记标准化第一潜变量被大数量纲属性独占。现象是 PLS 的潜变量个数选到 2 就收敛MAE 卡在 6 以上降不下去和论文的 3 到 4 差一截。原因是PLSRegression默认scaleTrue但如果你手动把 X 和 Y 拼在一起做了别的变换或者用的是自己实现的 PLS标准化这一步很容易漏。freq_sum这种属性的取值可能是几十而has_prefix只有 0 和 1不标准化时协方差完全由前者主导。解决方案是显式写scaleTrue并且在建模前打印每个属性的均值和标准差确认量纲已经被拉到同一水平。坑五混淆矩阵的标签顺序和聚类编号对不上。现象是准确率算出来只有 0.6 左右混淆矩阵的非对角元素散落但难度分的分布看起来完全正常。原因是 K-means 的簇编号是随机的第 0 簇不一定是难度最低的那一档直接拿簇编号和按分位切出来的档位比等于随机对齐。解决方式是先按每个簇内的平均难度分排序把簇编号重映射成 1 到 5再做混淆矩阵并且confusion_matrix里显式传labels[1,2,3,4,5]保证行列顺序固定。这五个坑里前三个属于数据层后两个属于模型层。我的排查顺序是先跑断言归一化、长度、日期连续性再看量纲最后才怀疑模型参数。大部分「复现不出来」的问题在数据层就能拦住。6. 二次开发把这套四问结构改成你自己的分析模板这份论文最值钱的地方不是那几个具体数值而是它给了一条可以整体搬走的结构清洗 → 单变量时序 → 多属性相关 → 多因变量回归 → 无监督分档 → 混淆矩阵验证。你换一个数据集只要满足「有时间维度、有对象属性、有比例型因变量」这三个条件这条链路基本可以原样复用。改造时我一般按下面这张表替换不动代码骨架只换数据接口原论文模块替换成你的场景需要改的地方num_reported日序列你的日活、订单量、访问量换列名确认是否有周效应9 个词属性你的对象特征商品、用户、内容重写features.py里的提取函数七档百分比你的比例型因变量分档转化率保证各行归一化加断言五档难度你的分层用户分层、内容分级重新跑肘部 轮廓系数定 k具体到一个最常见的改造——把「词」换成「商品」做难度分级。features.py里不再提取元音、重复字母改成提取价格区间、类目热度、标题长度、是否含促销词difficulty_score的权重向量从[1..7]换成你的转化档位权重PLS 的因变量从七档百分比换成各档购买转化率。骨架一行不用动preprocess.py里唯一要保留的是「归一化后加断言」和「时间序列补全」这两段这两段是通用的。# 换场景时只需要重写这个函数其余流程不变 def build_features(record: dict) - dict: title record[title] return { title_len: len(title), has_promo: int(any(w in title for w in [促销, 限时, 折扣])), price_log: np.log1p(record[price]), cat_heat: record[cat_heat], # 类目热度外部传入 }验证方法也很直接先跑一遍新数据的清洗断言确认没有非预期长度、没有归一化失败再画三张初检图看趋势是否符合业务直觉最后看混淆矩阵如果准确率低于 0.7大概率是属性提取漏了关键变量而不是模型选错了。我做完这套改造之后养成了一个习惯任何一份看起来「结构完整」的论文或模板我都先只跑清洗脚本把断言全打开跑通了再碰模型。因为数据层的静默错误比模型选错难查十倍而且它会一路污染到最后的结论里。希望这套拆法能帮到你少在数据那一步绕圈。本文还有配套的精品资源点击获取
简介这份2023年美国大学生数学建模竞赛C题的M奖Meritorious论文围绕Wordle游戏数据展开完整的建模与预测分析适合备战美赛的选手、需要参考英文论文写作规范的学生以及做数据分析课程设计或毕业设计的读者使用。全文为英文覆盖数据清洗与缺失值处理、9项单词属性的特征提取与Spearman相关性分析、时间序列区间预测、偏最小二乘回归PLS建模、K-means聚类结合肘部法则与轮廓系数确定词性与难度分组等环节并给出MAE在34、MSE在2224之间的误差评估与混淆矩阵热图还以单词EERIE为例算出4.635的难度得分并归入第五难度区间展示从聚类到评分的完整实现路径文末附有模型不足与不确定性讨论可作为完整赛题方案、论文结构与英文表达的学习范本。压缩包内共1个PDF文件约2.18MB为论文全文便于直接阅读与打印。目前已有170人学习下载适合希望快速理解M奖论文行文逻辑与建模思路的参赛者参考。1. 一份 25 页的 C 题论文卡住复现的往往不是公式去年有人拿着一份建模赛 C 题 M 奖论文来找我说只想把「3 月 1 日的结果数区间」这一问跑出来公式看得懂代码却怎么也复现不出[17363, 18994]这个区间。翻到数据那一节才发现卡点根本不在 ARIMA而在原始表里两个只有四个字母的单词和两处断崖式异常值上。这份论文做的是某每日五字母猜词游戏的公开数据用历史数据做时间序列区间预测、用九个词属性解释困难模式占比、用偏最小二乘预测某个词未来七档尝试次数的分布、再用 K-means 把单词按难度分级最后还留了一节讲数据里的有趣特征。它特别适合两类人一类是拿它当毕业设计或课程设计骨架、想找一条「清洗—特征—时序—回归—聚类—评估」完整链路的学生另一类是想看一份小样本多因变量建模怎么落地、又不想从零搭工程的开发者。下面按我拆它的顺序讲先看骨架再逐模块跑。2. 论文骨架拆解25 页目录怎么对应到一条可执行的数据流拿到这种论文我第一件事不是读模型是把目录抄成一张「章节 → 输入 → 输出 → 落地产物」的表。论文的章节编号本身已经暴露了作者的设计意图数据预处理单独成章放在最前面说明后四个问题共用同一份清洗结果第四问有趣特征放在模型之后说明它是对前面建模过程的观察沉淀不是独立任务。把这张表列出来你就知道哪一段该写脚本、哪一段只是结论陈述。论文章节对应任务输入输出II. Data preprocessing清洗 可视化359 行原始表12 列干净表、趋势图III. Problem 1时序预测 属性相关时间序列、9 个词属性3/1 预测区间、Spearman 系数IV. Problem 2多因变量回归词属性矩阵、七档百分比每档百分比的预测值、MAE/MSEV. Problem 3无监督难度分级词属性 预测百分比5 档难度、EERIE 得分、混淆矩阵VI. Problem 4特征观察全量数据若干规律性结论从数据流角度看整篇论文其实就一条主干一张宽表word, date, num_reported, num_hard, p1..p6, pX先被清洗成稳定序列再被拆成两个方向——时间维度走 ARIMA单词维度走属性提取 PLS最后两条线在第三问汇合用 PLS 预测出的七档百分比当权重算难度分再用 K-means 分档验证。项目如果要做成可复跑的代码仓我一般会按这个主干拆成四个文件wordle-c2023/ ├── data/ │ ├── raw/ # 原始 359 行数据只读不改 │ └── clean/ # 清洗后的中间结果可随时重建 ├── src/ │ ├── preprocess.py # 第二章纠错、异常值、百分比归一化 │ ├── forecast.py # 第三章ARIMA 定阶与区间预测 │ ├── features.py # 第三章9 个词属性 Spearman │ └── difficulty.py # 第四、五章PLS K-means 混淆矩阵 └── outputs/ ├── figures/ # 论文里的曲线图、饼图、散点图 └── tables/ # 预测区间、难度分、混淆矩阵这个结构看起来朴素但有两个好处raw目录只读保证任何时候都能从原始数据重跑一遍每一章对应一个脚本改模型不会牵一发动全身。论文里没有给出代码组织方式这是按常规工程习惯补的你完全可以把forecast.py和features.py合成一个只要保证清洗脚本先跑就行。2.1 第一刀清洗两个四字母单词和两处断崖式异常论文里写得很轻描淡写「发现有两个单词只有四个字母查证后修正为 clean 和 trash」。真跑起来这一步不能省因为后面所有基于「五字母」假设提取的属性首字母、尾字母、字母位置分布都会因为长度不对而错位。这类错误在公开数据集里很常见通常是录入时漏了字母处理方式就是人工核对后做一张映射表别用程序去猜。import pandas as pd import numpy as np df pd.read_csv(data/raw/wordle.csv) df[date] pd.to_datetime(df[date]) df df.sort_values(date).reset_index(dropTrue) # 1) 长度异常修正只对确认过的两个词做映射不做模糊匹配 typo_map {clen: clean, tash: trash} df[word] df[word].str.lower().str.strip().replace(typo_map) # 修正后必须全表满足五字母否则后面属性提取会静默出错 bad df[df[word].str.len() ! 5] assert bad.empty, f仍有非五字母单词{bad[word].tolist()} # 2) 断崖式异常用前后各一个点的均值替换 def fix_spike(frame, col, ts, window1): idx frame.index[frame[date] pd.Timestamp(ts)] if len(idx) 0: raise KeyError(f{ts} 不在数据集中) i idx[0] left, right frame.loc[i - window, col], frame.loc[i window, col] frame.loc[i, col] int(round((left right) / 2)) return frame df fix_spike(df, num_reported, 2022-11-30) # 2569 - (2373922628)/2 df fix_spike(df, num_hard, 2022-02-13) # 3249 - (931010343)/2这段代码有三个参数值得说清楚。typo_map必须是显式白名单我见过有人用「长度不等于 5 就丢弃」的写法结果把两个本来能救回来的样本扔了359 行变成 357 行后面的时序长度就和论文对不上。window1表示只取左右各一个邻居因为这两处异常都是单点尖刺窗口开大会把正常波动也抹平。round之后再取整是因为论文里的 23184 和 9827 都是四舍五入的结果不取整会在后续画图时出现小数刻度和论文图对不上。异常值为什么不能直接删因为这是时间序列删掉一行等于把 11 月 30 日和 12 月 1 日强行接在一起ARIMA 的滞后项会因此错位一天。均值替换虽然粗糙但保住了时间轴的连续性这是这套论文在数据量只有 359 行时比较务实的选择。2.2 百分比归一化98% 和 102% 都得变成 100%第二个坑更隐蔽。论文提到七档百分比之和在数据集中出现了 98%、99%、100%、101%、102% 五种情况。如果你直接拿原始百分比做 PLS 的因变量模型会学到一种假信号某些日期的「总概率质量」偏大或偏小但这个偏差其实来自四舍五入和平台统计口径不是玩家行为差异。pct_cols [p1, p2, p3, p4, p5, p6, pX] # 每行七档百分比求和再逐列除以总和强制归一 row_sum df[pct_cols].sum(axis1) df[pct_cols] df[pct_cols].div(row_sum, axis0) * 100 # 归一后校验每行之和与 100 的偏差应小于 1e-6 err (df[pct_cols].sum(axis1) - 100).abs().max() assert err 1e-6, f归一化未收敛最大偏差 {err:.6f} # 困难模式占比用人数直接相除比用百分比相减更稳 df[hard_ratio] df[num_hard] / df[num_reported] * 100这里div(row_sum, axis0)是按行广播axis0表示除数是每一行的标量容易写反成axis1写反了就是按列归一结果是每一列加起来等于 100看起来也「像对的」但完全破坏了样本间的比例关系属于典型的静默错误。断言里的1e-6是浮点容差太小会因为二进制表示误差误报太大又查不出问题。困难模式占比我没有用p1..pX那七列去推而是直接用num_hard / num_reported原因是论文第一问要分析的就是「困难模式得分百分比」和词属性的关系用人数相除得到的比例不受归一化舍入影响。这一步做完后面第三问的 Spearman 相关才有意义。2.3 可视化初检三张图其实是在给模型做体检论文第二章给了三组图七档百分比堆叠图、未成功挑战人数散点图、困难模式通关占比趋势图。很多人把这三张图当装饰其实它们各自对应一个建模假设的检验。堆叠图里黄色区域第四次尝试成功占比最大说明分布的重心在 3 到 5 次之间这直接决定了后面 PLS 的因变量不能当成七个独立回归硬拟合而要接受它们之间存在此消彼长的约束。未成功散点图里绝大多数日期低于 5%只有几十天超过 10%说明数据里存在少量「难词日」这些点在 K-means 分档时会被推到高档位如果直接当噪声剔掉难度分档的边界会整体左移。困难模式占比趋势图是单调上升的这一条最关键它意味着hard_ratio不是一个平稳序列如果你把它当因变量做时序回归必须先去趋势否则模型会把「时间」误当成「词属性」的解释力。我一般会在清洗脚本最后加一段快速出图只用来肉眼确认不追求论文里的排版import matplotlib.pyplot as plt fig, axes plt.subplots(1, 3, figsize(15, 4)) axes[0].stackplot(df[date], [df[c] for c in pct_cols]) axes[0].set_title(七档尝试次数占比归一后) axes[1].scatter(df[date], df[pX], s8) axes[1].axhline(5, ls--, cr) # 5% 参考线 axes[1].set_title(未成功挑战占比) axes[2].plot(df[date], df[hard_ratio]) axes[2].set_title(困难模式占比趋势) plt.tight_layout() plt.savefig(outputs/figures/preprocess_check.png, dpi150)axhline(5)那条红线不是论文原图里的是我自己加的参考线用来快速判断有多少天越过了 5% 这个量级。三张图放在一张画布上跑完清洗脚本顺手看一眼比等到建模阶段报错再回头查要省时间得多。3. 建模主干ARIMA 区间预测与 PLS 多因变量回归怎么落地清洗完之后论文的建模部分其实分成两条互不干扰的线。一条是时间维度只关心num_reported这个单变量序列目标是给出 2023 年 3 月 1 日的区间另一条是单词维度把每个词压成九个属性去解释困难模式占比再去预测七档百分比。两条线在代码上应该完全解耦因为它们的样本单位不一样前者是「天」后者是「词」混在一个脚本里最容易出现的错误就是把日期索引和单词索引对错行。3.1 ARIMA 定阶为什么最后落在低阶 (p,d,q)论文用 ARIMA(p,d,q) 做区间预测给出的结果是 [17363, 18994]。这三个参数里d 是差分阶数作用是让序列平稳p 是自回归项数q 是滑动平均项数。论文正文推导了(1-ΣφᵢLⁱ)(1-L)^d Xₜ (1ΣθᵢLⁱ)εₜ这个表达式落到代码里其实就是一次网格搜索加 AIC 选优。from statsmodels.tsa.arima.model import ARIMA # 以日期为索引缺失日期补全后再插值避免时间轴断裂 y df.set_index(date)[num_reported].asfreq(D).interpolate() best_aic, best_order, best_model np.inf, None, None for p in range(0, 4): for d in range(0, 3): for q in range(0, 4): try: m ARIMA(y, order(p, d, q), enforce_stationarityFalse, enforce_invertibilityFalse).fit() except Exception: continue # 不收敛的组合直接跳过 if m.aic best_aic: best_aic, best_order, best_model m.aic, (p, d, q), m print(最优阶数, best_order, AIC, round(best_aic, 2))参数上enforce_stationarityFalse和enforce_invertibilityFalse是为了让搜索空间更宽代价是可能选出物理意义不强的参数组合所以实战里我会把搜索上限压在 p4、q4。asfreq(D)加interpolate()这两步不能省因为 359 天里存在平台未记录的空档索引不连续时 statsmodels 的滞后计算会按行号而不是按日期走预测出来的区间会偏窄。拿到最优模型后预测区间这样出forecast best_model.get_forecast(steps(pd.Timestamp(2023-03-01) - y.index[-1]).days) pred forecast.predicted_mean.loc[2023-03-01] ci forecast.conf_int(alpha0.05).loc[2023-03-01] # 95% 置信区间 print(f点预测 {pred:.0f}区间 [{ci.iloc[0]:.0f}, {ci.iloc[1]:.0f}])alpha0.05对应 95% 置信水平论文里的 [17363, 18994] 就是这个口径下的结果。这里有个容易混淆的点置信区间和预测区间不是一回事前者描述的是「均值的估计误差」后者还要叠加单点观测的随机波动宽度更大。论文在正文里两个都建了代码上分别对应get_forecast().conf_int()和get_prediction().conf_int()如果你复现出来的区间比论文窄一截先检查是不是只取了均值置信区间。3.2 九个词属性怎么挑Spearman 比 Pearson 更适合这个场景论文提取了 9 个独立变量对困难模式得分百分比做 Spearman 相关分析结论是「这些属性有一定影响」。九个属性论文没有全部列出名称但从上下文能推断出大致方向词长这里恒为 5可作占位、重复字母个数、元音个数、字母频率总和、首字母频率、尾字母频率、是否含常见前缀、是否含常见后缀、字母分布熵。这套属性在猜词类游戏分析里是常规做法因为玩家猜词的难度主要来自「字母是否常见」和「结构是否有提示」。为什么用 Spearman 而不是 Pearson因为七档百分比是有界的且分布明显右偏第四次尝试那一档占比最大Pearson 假设线性关系和近似正态在这种有界比例数据上容易被极端值带偏。Spearman 只看秩对单调非线性更稳。from scipy.stats import spearmanr feat_cols [dup_cnt, vowel_cnt, freq_sum, first_freq, last_freq, has_prefix, has_suffix, entropy, uniq_cnt] rows [] for col in feat_cols: rho, pval spearmanr(df[col], df[hard_ratio]) rows.append({attr: col, rho: round(rho, 3), p: round(pval, 4)}) corr pd.DataFrame(rows).sort_values(rho, keyabs, ascendingFalse) print(corr)rho的绝对值超过 0.3 基本可以认为有中等程度关联p小于 0.05 说明这个关联不太可能是偶然。这里要注意九个属性彼此之间可能高度相关比如vowel_cnt和entropy单变量相关分析只能筛出候选不能直接下「因果」结论。论文的表述是「some effect」这个措辞是准确的因为后面 PLS 用的是全部九个属性一起进模型靠潜变量去吸收共线性而不是靠单变量筛选。我实际跑这类数据时会在相关分析后加一步 VIF 检查如果某个属性 VIF 超过 10就在 PLS 里降低它的权重或者干脆剔除否则潜变量个数会虚高。3.3 PLS 回归样本 359 行、七个因变量普通最小二乘为什么不够用第二问要把七个百分比1 到 6 次成功外加 X 未成功同时作为因变量预测自变量是九个词属性。这是一个典型的多因变量小样本问题样本量 359自变量 9 个因变量 7 个而且七个因变量之和恒等于 100存在完全共线性——你不可能同时独立地预测七个数因为它们被一个约束绑死了。普通多元线性回归在这里会直接翻车系数矩阵不满秩。论文选 PLS偏最小二乘是对的。PLS 的思路是先对 X 和 Y 同时做潜变量提取类似 PCA但提取方向要最大化与 Y 的协方差再在潜变量空间里做回归。它把主成分分析、多元线性回归和误差分析揉在一起对样本量适中、变量间共线性强的数据特别友好。from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import cross_val_predict, KFold from sklearn.metrics import mean_absolute_error, mean_squared_error X df[feat_cols].values Y df[pct_cols].values # 5 折交叉验证选潜变量个数避免用测试集调参 kf KFold(n_splits5, shuffleTrue, random_state42) best_n, best_mae 1, np.inf for n in range(1, 8): pls PLSRegression(n_componentsn, scaleTrue) Y_pred cross_val_predict(pls, X, Y, cvkf) mae mean_absolute_error(Y, Y_pred) if mae best_mae: best_n, best_mae n, mae final PLSRegression(n_componentsbest_n, scaleTrue).fit(X, Y) Y_hat final.predict(X) print(潜变量个数, best_n, MAE, round(mean_absolute_error(Y, Y_hat), 2), MSE, round(mean_squared_error(Y, Y_hat), 2))scaleTrue是关键参数它会对 X 和 Y 都做标准化。九个属性的量纲差别很大——freq_sum可能是几十has_prefix只有 0/1——不标准化的话量纲大的属性会独占第一个潜变量。n_components论文没有明说取了几个但常规做法是在 1 到 5 之间用交叉验证选选太大就会过拟合交叉验证 MAE 会先降后升。论文给出的 MAE 在 3 到 4 之间、MSE 在 22 到 24 之间这个量级是合理的七档百分比的均值大概在 14 左右MAE 3 到 4 意味着平均每档误差三到四个百分点。要注意 MAE 和 MSE 的单位是「百分点」不是比例所以 22 到 24 的 MSE 开根号后 RMSE 大约 4.7和 MAE 同量级说明误差分布没有特别重的长尾。4. 难度分级与验证K-means 分档、加权打分和 0.941 的准确率第三问是整个论文里最容易写得好看、也最容易复现翻车的一问。任务本身不复杂单词没有现成的难度标签所以用无监督学习分档分完档之后用 PLS 预测的七档百分比算一个加权难度分最后用混淆矩阵验证分档的准确率。三个环节环环相扣任何一环的标签顺序错了最后那个 0.941 就出不来。4.1 聚类数怎么定肘部法则看拐点轮廓系数看分离度K-means 最大的问题是 k 要自己给。论文用了两个指标交叉验证肘部法则和轮廓系数。肘部法则看的是簇内平方和inertia随 k 变化的曲线拐点处就是「再加一个簇带来的收益开始变小」的位置轮廓系数衡量的是每个点到本簇的紧密度和到最近邻簇的分离度取值 -1 到 1越高越好。from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler from sklearn.metrics import silhouette_score Xw StandardScaler().fit_transform(df[feat_cols]) # 聚类前必须标准化 inertia, sil [], {} for k in range(2, 11): km KMeans(n_clustersk, n_init20, random_state42).fit(Xw) inertia.append(km.inertia_) sil[k] silhouette_score(Xw, km.labels_) # 肘部相邻两点斜率下降最剧烈的位置 drops np.diff(inertia) / inertia[:-1] elbow_k int(np.argmin(drops)) 2 # 2 是因为 k 从 2 开始 print(肘部候选 k , elbow_k, 轮廓系数 , round(sil[elbow_k], 3))n_init20表示用 20 组不同初始质心各跑一次取最优K-means 对初始化敏感默认的 10 次在小样本上偶尔会掉进局部最优。random_state42是为了可复现换了随机种子簇编号会变但分档结构一般不变——这一点很重要因为后面混淆矩阵的标签对齐依赖簇编号固定种子能省掉一堆麻烦。论文最后定的是五档难度这个数字通常不是肘部法则直接给出的而是肘部候选常见是 3 到 4和轮廓系数峰值之间折中再结合「五档在业务上好解释」这个约束定下来的。我一般会把 inertia 和轮廓系数的曲线都画出来放在一起看两个指标都支持的那个 k 才靠谱。4.2 难度分怎么算EERIE 的 4.635 是加权期望分完档只是有了「组」组内怎么排序、怎么给出一个连续分数论文用的是加权公式把 PLS 预测出的七档百分比当作权重乘上对应的尝试次数求和得到一个期望尝试次数值越大表示越难。# p1..p6 对应 1~6 次成功pX 是未成功权重记为 7 weights np.array([1, 2, 3, 4, 5, 6, 7]) def difficulty_score(pct_row): pct np.asarray(pct_row, dtypefloat) pct pct / pct.sum() # 防止归一化残差 return float((pct * weights).sum()) eer df.loc[df[word] eerie, pct_cols].values[0] print(EERIE 难度分 , round(difficulty_score(eer), 3)) # 论文中为 4.635权重 7 是给未成功那一档的这个取值不是唯一的取 7 相当于假设「失败比第六次成功再难一档」你也可以取 6.5 或 8难度分的绝对值会变但排序基本稳定。EERIE 这个词有重复字母、元音密集属于典型的「看着简单猜着难」4.635 落在五档里的高档位和论文结论一致。这一步有个容易忽略的细节PLS 预测出来的百分比可能为负或者和不为 100所以pct / pct.sum()这行归一化不能省否则加权和会出现负数或者超过 7 的异常值。4.3 混淆矩阵 0.941验证的不是聚类是「预测能不能落回原档」论文用混淆矩阵热图评估模型准确率 0.941。这里要理解清楚它在验证什么K-means 本身是无监督的没有真标签所以混淆矩阵的标签来源只能是「用预测百分比重新算难度分再按难度区间反推档位」和原始聚类结果对比。换句话说它验证的是「PLS 预测 → 加权打分 → 分档」这条链路和直接聚类的结果是否一致一致性高说明 PLS 的预测没有把词的结构信息丢掉。from sklearn.metrics import confusion_matrix, accuracy_score # 用全量数据重新算难度分并按分位切五档 scores np.array([difficulty_score(row) for row in df[pct_cols].values]) bins np.quantile(scores, [0, 0.2, 0.4, 0.6, 0.8, 1.0]) pred_level np.digitize(scores, bins[1:-1]) 1 # true_level 来自 K-means 的簇编号按簇内均分排序后映射到 1~5 cm confusion_matrix(true_level, pred_level, labels[1, 2, 3, 4, 5]) acc accuracy_score(true_level, pred_level) print(cm) print(准确率, round(acc, 3)) # 论文中为 0.941np.quantile用等频分箱保证五档样本量接近这比等宽分箱更稳因为难度分是右偏的等宽分箱会把大部分词塞进低档。labels[1,2,3,4,5]显式指定标签顺序防止混淆矩阵行列顺序错位——这是复现这类结果最常见的翻车点代码跑通但数值对不上八成是标签顺序反了。准确率 0.941 意味着 359 个词里大约 21 个被分到了相邻档位这个误差水平在无监督分档里算相当好。看混淆矩阵时重点看对角线附近如果误差主要集中在 2 和 3 之间说明这两档的边界本来就模糊属于正常如果误差散落在对角线上那就要怀疑属性提取有问题。4.4 第四问的特征观察几个值得单独记下来的规律第四问没有建模是对全量数据的观察总结论文列了四条时间上的相关性、首字母频率、尝试次数的近似正态分布、以及某些特殊词缀的异常。这几条其实都能从前面清洗后的数据里直接验证。时间相关性指的是结果数序列存在周内波动周末参与人数偏高这也是 ARIMA 里 d 取小的原因之一过度差分会把周效应也差掉。首字母频率说的是以常见字母开头的词玩家第一次猜中的概率更高这会直接影响p1这一档的百分比分布。尝试次数近似正态这一点和第二章的堆叠图对得上重心在第四次。特殊词缀异常是指含常见后缀的词反而更难因为玩家容易先入为主地往错误方向猜。这四条结论的价值不在于多深奥而在于它们都是「建模过程中顺手能看到」的副产物。我拆这类论文时习惯把第四问当成检查清单如果你复现出来的数据也能观察到这四条说明前面的清洗和特征提取方向是对的如果一条都看不到先回头查归一化和异常值处理。5. 避坑与排查复现这套代码时最容易翻车的五处这一章是我自己踩过的坑按「现象 → 原因 → 解决」写。这五个问题有一个共同特点代码不报错结果看起来正常但数值就是和论文对不上。5.1 数据层的三个坑坑一百分比归一化写成了按列归一。现象是模型能跑MAE 也不大但七档百分比的预测值加起来不等于 100强行解释成「模型误差」。原因是df[pct_cols].div(row_sum, axis0)里的axis写反了写成了axis1变成每个列向量除以一个标量序列实际效果是按列缩放。解决方式很简单在归一化后加一行断言(df[pct_cols].sum(axis1) - 100).abs().max() 1e-6让它在清洗阶段就炸出来而不是等到建模阶段猜。坑二异常值替换时用错了邻居。现象是替换后序列在异常点附近出现一个明显的台阶ARIMA 的预测区间整体偏移。原因是代码里用了df.loc[i-1]和df.loc[i1]取邻居但这个 DataFrame 的索引在排序后没有重置loc取的是标签不是位置取到了别的日期。解决方式是清洗开始时统一reset_index(dropTrue)替换函数里改用iloc[i-window]和iloc[iwindow]并且把替换前后的值打印出来核对2569 变 23184、3249 变 9827这两个数对上了再往下走。坑三时间序列没补全就进 ARIMA。现象是模型能拟合但预测到 3 月 1 日的步数和论文差几天区间宽度明显偏窄。原因是原始数据里有些日期平台没记录直接set_index(date)后索引不连续statsmodels 按行号推算滞后把不连续的日子当成了连续。解决方式是asfreq(D)补齐日历日再interpolate()插值插值方法用线性就够因为缺的是零星单日不是整段缺失。5.2 模型层的两个坑坑四PLS 忘记标准化第一潜变量被大数量纲属性独占。现象是 PLS 的潜变量个数选到 2 就收敛MAE 卡在 6 以上降不下去和论文的 3 到 4 差一截。原因是PLSRegression默认scaleTrue但如果你手动把 X 和 Y 拼在一起做了别的变换或者用的是自己实现的 PLS标准化这一步很容易漏。freq_sum这种属性的取值可能是几十而has_prefix只有 0 和 1不标准化时协方差完全由前者主导。解决方案是显式写scaleTrue并且在建模前打印每个属性的均值和标准差确认量纲已经被拉到同一水平。坑五混淆矩阵的标签顺序和聚类编号对不上。现象是准确率算出来只有 0.6 左右混淆矩阵的非对角元素散落但难度分的分布看起来完全正常。原因是 K-means 的簇编号是随机的第 0 簇不一定是难度最低的那一档直接拿簇编号和按分位切出来的档位比等于随机对齐。解决方式是先按每个簇内的平均难度分排序把簇编号重映射成 1 到 5再做混淆矩阵并且confusion_matrix里显式传labels[1,2,3,4,5]保证行列顺序固定。这五个坑里前三个属于数据层后两个属于模型层。我的排查顺序是先跑断言归一化、长度、日期连续性再看量纲最后才怀疑模型参数。大部分「复现不出来」的问题在数据层就能拦住。6. 二次开发把这套四问结构改成你自己的分析模板这份论文最值钱的地方不是那几个具体数值而是它给了一条可以整体搬走的结构清洗 → 单变量时序 → 多属性相关 → 多因变量回归 → 无监督分档 → 混淆矩阵验证。你换一个数据集只要满足「有时间维度、有对象属性、有比例型因变量」这三个条件这条链路基本可以原样复用。改造时我一般按下面这张表替换不动代码骨架只换数据接口原论文模块替换成你的场景需要改的地方num_reported日序列你的日活、订单量、访问量换列名确认是否有周效应9 个词属性你的对象特征商品、用户、内容重写features.py里的提取函数七档百分比你的比例型因变量分档转化率保证各行归一化加断言五档难度你的分层用户分层、内容分级重新跑肘部 轮廓系数定 k具体到一个最常见的改造——把「词」换成「商品」做难度分级。features.py里不再提取元音、重复字母改成提取价格区间、类目热度、标题长度、是否含促销词difficulty_score的权重向量从[1..7]换成你的转化档位权重PLS 的因变量从七档百分比换成各档购买转化率。骨架一行不用动preprocess.py里唯一要保留的是「归一化后加断言」和「时间序列补全」这两段这两段是通用的。# 换场景时只需要重写这个函数其余流程不变 def build_features(record: dict) - dict: title record[title] return { title_len: len(title), has_promo: int(any(w in title for w in [促销, 限时, 折扣])), price_log: np.log1p(record[price]), cat_heat: record[cat_heat], # 类目热度外部传入 }验证方法也很直接先跑一遍新数据的清洗断言确认没有非预期长度、没有归一化失败再画三张初检图看趋势是否符合业务直觉最后看混淆矩阵如果准确率低于 0.7大概率是属性提取漏了关键变量而不是模型选错了。我做完这套改造之后养成了一个习惯任何一份看起来「结构完整」的论文或模板我都先只跑清洗脚本把断言全打开跑通了再碰模型。因为数据层的静默错误比模型选错难查十倍而且它会一路污染到最后的结论里。希望这套拆法能帮到你少在数据那一步绕圈。本文还有配套的精品资源点击获取
