简介：这份2023年美赛C题M奖论文以Wordle游戏数据为研究对象，面向美赛参赛者、数据分析学习者及需要毕业设计参考的读者，提供一套从数据预处理到建模评估的完整研究范例。全文围绕时间序列预测、不同尝试次数得分百分比预测、单词词性与难度分类三个问题展开，涉及缺失值与异常文本处理、9个单词属性提取及Spearman相关性分析。模型部分采用时间序列区间预测、偏最小二乘回归与K-means聚类，并结合肘部法则、轮廓系数确定分类数，给出MAE约3至4、MSE约22至24及准确率0.941等评估结果。文中还讨论了数据质量、模型假设等不确定性因素，并给出后续研究展望，便于读者借鉴论文组织方式与论证逻辑。资源包为1个PDF文件，大小2.18MB，便于直接阅读与引用；已有170人学习下载，适合用作赛题复盘、论文结构模仿和建模方法查阅。

1. 从一份 359 行的 Wordle 表开始：它到底能替你做掉哪些建模活

做美赛 C 题那年，很多人卡在同一个地方：数据给到手了，四问也读懂了，可真要落笔写模型时，发现时间序列、多因变量回归、无监督聚类三套东西得在一篇论文里串起来，串不好就是四段互不相干的作业拼盘。这份 2023 年美赛 C 题的 M 奖论文正好是一份完整拆解样本：359 行有效数据、12 列字段，从异常值替换、百分比归一化，一路做到 ARIMA 区间预测、9 个词属性的 Spearman 相关、PLS 多因变量回归，最后用 K-means 把单词难度分成五档，给出 EERIE 4.635 分的结论。它解决的不是「Wordle 怎么玩」，而是「一份带时间维度、多因变量、无标签的数据，怎么在一篇 25 页论文里把四问闭环」。适合正在做美赛 C 题、或者拿这类多方法串联题目当毕业设计的人。

2. 数据清洗：把 12 列原始表洗成能进模型的矩阵

这份数据的第一眼印象是「干净得可疑」——359 行、12 列，简介里明说没有缺失值。但真跑起来，四个地方会连续翻车：两个只有四个字母的单词、五条加不到 100% 的百分比行、两个把曲线戳出尖刺的异常点、以及一个必须补成等间隔的日期索引。清洗没做干净，后面 ARIMA 的残差检验一定过不了。

2.1 先看三件事：词长、百分比之和、日期索引

拿到表别急着画图，先做三个断言。词长不是 5 的行，说明录入时被截断了；七列百分比之和不等于 1 的行，说明原始数据做过四舍五入；日期列不连续，说明有缺日或者重复日，直接喂给时间序列模型会报错或者悄悄错位。

论文里那两个词是clen和tash，翻游戏记录后分别修正为clean和trash。这个细节很值得学——它不是随手改的，是拿「Wordle 每期答案必须是一个真实的五字母英文单词」这条硬约束反推出来的。做毕业设计时如果数据源是你自己爬的，同样得先立一条这样的硬约束，再谈清洗。

import pandas as pd import numpy as np # 359 行 × 12 列，列名按常见结构命名 df = pd.read_csv("wordle.csv") print(df.shape) # (359, 12) # 断言 1：词长必须为 5 bad_word = df[df["Word"].str.len() != 5] print(bad_word[["Date", "Word"]]) # 会打出 clen / tash 两行 fix_map = {"clen": "clean", "tash": "trash"} df["Word"] = df["Word"].replace(fix_map) # 断言 2：七列百分比之和是否都等于 1 pct_cols = ["1 try", "2 tries", "3 tries", "4 tries", "5 tries", "6 tries", "7 or more tries (X)"] s = df[pct_cols].sum(axis=1) print(s.round(2).value_counts()) # 0.98 / 0.99 / 1.00 / 1.01 / 1.02 五种 # 断言 3：日期是否等间隔 df["Date"] = pd.to_datetime(df["Date"]) gap = df["Date"].diff().dt.days.value_counts() print(gap) # 正常应该只有 1

这段代码的价值不在语法，在于它把「肉眼检查」变成了「可重复的断言」。参数上唯一要留意的是pct_cols的列名——不同年份、不同来源的数据表列名可能是1 try、1_try、Score1，进模型前必须统一，否则后面 PLS 的Y矩阵会静默少列。日期那一步如果value_counts()出现了 2 或者 3，说明有跳日，得先用asfreq("D")补齐再插值，不能直接set_index了事。

2.2 异常值不是删掉，是用两侧均值顶掉

论文里画了两条曲线：Number of reported results和Number in hard mode随时间的变化。两条曲线上各有一个明显的尖刺——2022/11/30 的答题人数只有 2569，前后两天分别是 23739 和 22628；2022/2/13 的困难模式人数只有 3249，前后是 9310 和 10343。这两个点如果留着，ARIMA 的方差会被抬高，预测区间会被撑得很宽，最后给出的区间就没参考价值了。

处理方式很朴素：取前后两点的均值，四舍五入。(23739 + 22628) / 2 = 23183.5，取 23184；(9310 + 10343) / 2 = 9826.5，取 9827。这套做法在时间序列里叫「三点平滑」的最简形式，优点是保序、不引入新趋势，缺点是对连续多个异常点无能为力。359 行里只有两个孤立尖刺，用它足够了。

def smooth_by_neighbors(series, idx): """用前后两点的均值替换单点异常，idx 为行号""" left, right = series.iloc[idx - 1], series.iloc[idx + 1] return int(round((left + right) / 2)) # 定位异常点（也可以用 3σ 或者 IQR 自动找） mask = df["Number of reported results"] < 5000 for idx in df.index[mask]: df.loc[idx, "Number of reported results"] = smooth_by_neighbors( df["Number of reported results"], idx) mask = df["Number in hard mode"] < 5000 for idx in df.index[mask]: df.loc[idx, "Number in hard mode"] = smooth_by_neighbors( df["Number in hard mode"], idx)

参数说明：阈值 5000 不是拍脑袋定的，是因为正常日期的答题量在万级，2569 和 3249 明显掉了一个数量级。如果你的数据量级不同，换成「中位数 ± 3 倍 MAD」更稳：

med = df["Number of reported results"].median() mad = (df["Number of reported results"] - med).abs().median() mask = (df["Number of reported results"] - med).abs() > 3 * 1.4826 * mad

这一步的常见误用是把异常点直接drop掉。少一行看起来无所谓，但时间序列一旦断档，asfreq("D")会补出NaN，后面差分、ADF 检验全部受影响。替换比删除安全，这是我踩过一次才记住的。

2.3 百分比归一化：除法和补差值的取舍

七列百分比加起来有 98%、99%、100%、101%、102% 五种情况，源头是原始数据保留整数位后的舍入误差。不归一化，PLS 的Y矩阵每一行的总和不一致，回归出来的系数会带上一个跟单词无关的常数偏移。

论文用的是 Excel 里逐行按比例缩放。换成代码就是整行除以行和：

raw_sum = df[pct_cols].sum(axis=1) df[pct_cols] = df[pct_cols].div(raw_sum, axis=0) print(df[pct_cols].sum(axis=1).round(6).unique()) # 全部为 1.0

div(..., axis=0)的axis=0是关键——它表示按行广播，写错成axis=1会变成按列除，结果每一列加起来是 1，模型直接废掉。这个错我在第一次写的时候犯过，因为 pandas 的axis语义和 numpy 是反的，血泪经验。

还有一种写法是把误差全部补给X（七次以上未成功）那一列，理由是它本来就是残差项。两种写法在 100% 附近差别很小，但如果某一行是 98%，按比例缩放会把 1 try 到 6 tries 全部放大 2%，而补差值只动 X 一列。做难度分级时，X的占比本身就是难度信号，所以我一般会选按比例缩放，避免人为污染这个特征。

归一化之后再做一次分布检查：把 1 到 6 tries 画成堆叠面积图，正常形态应该是 3、4、5 次尝试占大头，X稳定在 5% 以下。论文里那张图就是这个结论，它还顺带说明了一件事——359 天里只有几十天的失败率超过 10%，说明官方对难度的控制相当稳。这个观察在后面 K-means 分档时可以直接拿来当合理性校验：如果聚出来的「最难档」对应的X占比只有 3%，那这个聚类结果八成有问题。

3. ARIMA 区间预测：从 ADF 检验到 3 月 1 日的 [17363, 18994]

第一问要的是「2023 年 3 月 1 日报告人数的预测区间」。注意是区间，不是点估计——这是很多人丢分的地方，只给一个数，等于没回答「vary daily」这件事。论文给的结果是 [17363, 18994]，宽度大约 1600，相对中位数约 9%，这个宽度是合理的，太窄说明方差被低估，太宽说明模型没抓住趋势。

3.1 为什么是 ARIMA，而不是线性回归或 LSTM

359 个观测点，日频，单变量，带明显趋势但没有强季节性——这个数据画像基本就是 ARIMA 的主场。线性回归得自己构造滞后特征，还得处理残差自相关；LSTM 在这个样本量下大概率过拟合，而且给不出解析形式的预测区间，只能靠 bootstrap 硬凑。

ARIMA(p, d, q) 里三个参数各有分工：p是自回归项数，刻画「昨天的人数影响今天」；d是差分阶数，把非平稳序列压成平稳；q是滑动平均项数，刻画历史冲击的残留影响。论文里用公式 (1) 到 (5) 把这套东西推了一遍，核心结论就一句：d=0时 ARIMA 退化成 ARMA，所以差分阶数是区分平稳与非平稳的开关。

选型上还有一个现实理由：ARIMA 的预测区间是模型内置的，get_forecast().conf_int()一行就能拿到，不用额外写代码。对于三天的比赛节奏，这个性价比很高。

3.2 三步定阶：ADF、差分、Ljung-Box

定阶不要凭感觉试 (1,1,1)、(2,1,2)，按顺序走三步，每一步都有明确的通过标准。

第一步，ADF 检验看原序列是否平稳。p 值大于 0.05 就拒绝「平稳」的原假设，需要差分。第二步，一阶差分后再做一次 ADF，通常到这一步就平稳了；如果还不行，d=2，但一般日频数据不需要。第三步，拟合之后对残差做 Ljung-Box 检验，p 值大于 0.05 说明残差已经没有自相关，模型把信息榨干了。

import pandas as pd from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.arima.model import ARIMA from statsmodels.stats.diagnostic import acorr_ljungbox y = df.set_index("Date")["Number of reported results"].asfreq("D") # 第一步：原序列 ADF print("raw p =", adfuller(y, autolag="AIC")[1]) # 第二步：一阶差分后 ADF print("diff p =", adfuller(y.diff().dropna(), autolag="AIC")[1]) # 第三步：拟合 + 残差白噪声检验 model = ARIMA(y, order=(1, 1, 1)).fit() print(model.summary()) print(acorr_ljungbox(model.resid, lags=[7, 14], return_df=True)) # 预测 2023-03-01 fc = model.get_forecast(steps=1) print("mean :", fc.predicted_mean.iloc[0]) print("95% CI:", fc.conf_int(alpha=0.05).iloc[0].tolist())

参数说明：autolag="AIC"让 statsmodels 自动选滞后阶数，比固定maxlag稳；lags=[7, 14]是因为日频数据一周、两周的自相关最值得看，只看 lag=1 会漏掉周内周期。如果 Ljung-Box 的 p 值小于 0.05，说明残差还有结构，优先调p和q，不要急着加d——差分过度会让序列过平稳，预测区间反而失真。

定阶的搜索空间不用开太大，p和q各取 0 到 3，一共 16 组，用 AIC 最小选一组就够。样本只有 359 个，高阶模型没有足够的数据支撑，这也是论文里没往深里做的原因。

3.3 置信区间和预测区间别混用

这两个概念在论文里都提到了，但很多复现的人只取了一个。区别在于：置信区间描述的是「均值的估计有多准」，预测区间描述的是「下一次实际观测会落在哪」。后者一定比前者宽，因为多了观测噪声这一项。

对于「3 月 1 日的报告人数」这个问题，要的是预测区间。statsmodels 的conf_int()默认给的是均值置信区间，想要预测区间得手动加方差：

import numpy as np from scipy import stats alpha = 0.05 mean = fc.predicted_mean.iloc[0] se_mean = (fc.conf_int(alpha=0.05).iloc[0, 1] - mean) / stats.norm.ppf(1 - alpha / 2) sigma = np.sqrt(model.params["sigma2"]) z = stats.norm.ppf(1 - alpha / 2) lower = mean - z * np.sqrt(se_mean ** 2 + sigma ** 2) upper = mean + z * np.sqrt(se_mean ** 2 + sigma ** 2) print(f"prediction interval: [{lower:.0f}, {upper:.0f}]")

sigma2是残差方差，se_mean是均值标准误，两者平方和开根号就是单次预测的标准差。跑出来的区间会比conf_int()宽一截，这才对得上论文里 [17363, 18994] 这个量级。

还有一个容易忽略的点：asfreq("D")之后如果中间有缺日，statsmodels 会带着NaN一起拟合，结果不报错但系数是错的。清洗阶段做的日期连续性断言，就是为这一步服务的。

4. 词属性到 PLS：9 个自变量怎么撑起 7 个百分比

第二问要预测一个未来单词在 (1,2,3,4,5,6,X) 七个类别上的百分比分布。这是一个「多输入多输出」问题：输入是单词的各种属性，输出是七个百分比。普通线性回归得训七个模型，还得处理七个输出之间的相关性——七个百分比加起来等于 1，天然负相关，分开训必然互相打架。

4.1 9 个词属性的提取口径

论文里提取了 9 个独立变量，并对它们和「困难模式得分百分比」做了 Spearman 相关分析，结论是这 9 个属性都有一定影响。具体是哪 9 个，原文没逐条列全，但按这个题的标准做法，通常是下面这组：

import numpy as np # freq: dict，字母 -> 在语料中出现的频率 def build_features(word, freq, pos_dist, word_rank): letters = list(word) return { "dup_letters": len(letters) - len(set(letters)), "vowels": sum(c in "aeiou" for c in letters), "rare_cnt": sum(freq.get(c, 0) < 0.02 for c in letters), "mean_freq": float(np.mean([freq.get(c, 0) for c in letters])), "first_freq": freq.get(letters[0], 0), "last_freq": freq.get(letters[-1], 0), "pos_entropy": sum( -sum(p * np.log(p + 1e-9) for p in pos_dist[i].values()) for i in range(5)), "double": int(any(letters[i] == letters[i + 1] for i in range(4))), "log_rank": np.log1p(word_rank.get(word, 10000)), } X = pd.DataFrame([build_features(w, freq, pos_dist, rank) for w in df["Word"]]) print(X.shape) # (359, 9)

参数上最需要注意的是pos_dist——每个位置上字母的出现分布，必须用跟当前数据集同期的语料统计，不能用今天的语料。否则pos_entropy会带上时间泄漏，模型在训练集上好看，一预测未来就崩。这是做毕业设计时最隐蔽的一个坑。

4.2 Spearman 先筛一遍再进 PLS

九个属性之间不是独立的：mean_freq和first_freq高度相关，rare_cnt和mean_freq是反向关系。直接进回归会有共线性，系数符号可能反直觉。

论文用 Spearman 而不是 Pearson，是因为这些属性和百分比之间的关系未必线性，而且百分比里可能有离群值。Spearman 只看秩，对异常值不敏感。做法是算每个属性和目标（困难模式百分比）的秩相关系数和 p 值，把 p 值大于 0.05 的剔除。

from scipy.stats import spearmanr target = df["% in hard mode"] rows = [] for col in X.columns: rho, p = spearmanr(X[col], target) rows.append({"attr": col, "rho": round(rho, 3), "p": round(p, 4)}) corr_table = pd.DataFrame(rows).sort_values("rho", key=abs, ascending=False) print(corr_table)

结果表里如果某个属性的rho接近 0 且p > 0.05，说明它和目标没关系，删掉。但注意别一次删太多——九个属性删到只剩三个，模型解释性会变差，论文里也不好写。我的做法是保留|rho| > 0.1或者p < 0.1的，宁可留一点冗余，PLS 本身对共线性有容忍度。

4.3 n_components 怎么定，MAE 3~4 才算过关

PLS 的核心是把 X 和 Y 同时投影到低维空间，找一组潜变量，让投影后的协方差最大。它相当于「PCA + 多元线性回归」的合体，特别适合样本量中等、自变量有共线性的场景——359 个样本，9 个输入，7 个输出，正好落在它的舒适区。

from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import cross_val_predict from sklearn.metrics import mean_absolute_error, mean_squared_error Y = df[pct_cols].values * 100 # 转成百分数，量纲对齐论文 best = None for k in range(1, 7): pls = PLSRegression(n_components=k, scale=True) pred = cross_val_predict(pls, X, Y, cv=5) mae = mean_absolute_error(Y, pred) mse = mean_squared_error(Y, pred) print(f"k={k} MAE={mae:.3f} MSE={mse:.3f}") if best is None or mse < best[1]: best = (k, mse) k, _ = best pls = PLSRegression(n_components=k, scale=True).fit(X, Y) print("chosen k =", k, " coef shape =", pls.coef_.shape)

scale=True必须开。九个属性量纲完全不同——dup_letters是 0 到 4 的整数，pos_entropy可能是 5 到 10 的浮点，不标准化的话 PLS 会把权重大头给数值大的特征。

n_components控制潜变量个数，太小欠拟合，太大会把噪声也学进去。用 5 折交叉验证的 MSE 最小来选，比看训练集误差靠谱得多。论文里给出的验收标准是 MAE 落在 3 到 4 之间、MSE 落在 22 到 24 之间，单位是百分点。也就是说，七个百分比里每一项的平均绝对误差不超过 4 个百分点。这个精度对于「预测一个还没出现的单词的答题分布」已经够用了——毕竟人类自己猜一个陌生单词的难度，误差也差不多这个量级。

如果跑出来 MAE 超过 6，先查三件事：百分比有没有归一化、Y有没有乘 100、特征里有没有时间泄漏。这三个是最常见的翻车点，跟模型本身没关系。

5. 五个必踩的坑：从 98% 的百分比之和到 K-means 不缩放

这一节是我把论文流程完整复现一遍之后攒下来的踩坑记录，每条都按「现象 → 原因 → 解决」写。有些坑原文一笔带过，但真跑代码的时候会卡住半天。

5.1 七列百分比加起来是 98%，PLS 系数全偏

现象：模型训完，七个输出的残差均值不接近 0，而是统一偏了一个正数或负数，MAE 稳定在 5 以上降不下去。

原因：原始数据里 359 行中有相当一部分百分比之和是 98%、99%、101%、102%，没有归一化就乘 100 进了Y矩阵。PLS 拟合的是绝对数值，行和不为 100 相当于给每一行加了一个跟单词无关的偏移，这个偏移被模型当成了信号的一部分。

解决：进模型前强制df[pct_cols].div(df[pct_cols].sum(axis=1), axis=0)，归一化后再乘 100。归一化之后还要再断言一次sum == 1.0，浮点误差范围内允许 1e-6 的偏差。

5.2 ADF 检验 p 值 0.9，差分两次后区间宽得没法用

现象：原序列 ADF 的 p 值接近 1，一阶差分后降到 0.03，觉得还不够稳，又做了一次差分，结果预测区间宽度翻倍。

原因：一阶差分已经让序列平稳了，第二次差分属于过度差分。过度差分会在序列里引入人为的负自相关，ARIMA 为了拟合这个假结构，会把方差估计抬高，预测区间自然被撑宽。

解决：ADF 的判定标准是 p < 0.05 就停，不要追求更小的 p 值。差分之后再看一眼 ACF 图，如果 lag=1 的自相关是明显负数，基本可以确定差分过头了。d的上限设成 2，实际绝大多数日频数据一次就够。

5.3 日期没补全，ARIMA 悄悄算错

现象：模型能跑通，summary()也没有警告，但预测出来的数跟实际差了 20% 以上。

原因：原始表里的日期有跳日，set_index之后没有asfreq("D")，statsmodels 会按行索引当成等间隔处理。中间缺的那一天被吞掉，所有滞后项的位置都往前挪了一位，模型学到的是错位的自相关。

解决：set_index("Date").asfreq("D")，然后检查isna().sum()。如果补出了缺失值，用线性插值填，别用前向填充——前向填充会让连续几天的值一模一样，差分会得到一串 0，同样污染模型。

5.4 K-means 不标准化，聚类结果被一个特征绑架

现象：聚出来的五档难度，每一档的log_rank高度集中，其他属性在档内几乎没差异，但业务上看这个分档毫无意义。

原因：九个属性里pos_entropy的取值范围可能是 0 到 10，double只有 0 和 1，欧氏距离会被数值大的特征主导。K-means 的本质是最小化点到质心的欧氏距离，量纲不统一，等于只拿一个特征在聚类。

解决：StandardScaler先做 z-score，再进 KMeans。标准化之后每个特征的方差都是 1，距离才有可比性。这一步跟 PLS 里的scale=True是同一个道理，但很多人做聚类时会忘——因为聚类看起来「无参数」，容易掉以轻心。

5.5 用训练集准确率 0.941 下结论，被追问就答不上来

现象：论文里写分类准确率 0.941，答辩或者评审一问「这个准确率是在哪个集合上算的」，答不出来。

原因：K-means 是无监督的，本身没有标签，所谓「准确率」一定是拿聚类结果跟某个人工标注的难度档对比出来的。如果这个对比用的是全量数据，而且 K 值也是在全量数据上选的，那就存在信息泄漏。

解决：K 值选择（肘部法则、轮廓系数）和最终评估要分开。常见做法是先把数据按 8:2 切分，在训练集上选 K、训质心，在测试集上算混淆矩阵和准确率。如果样本量太小不方便切，至少要在论文里写清楚「准确率的计算口径」，别只甩一个数字。

6. 轮廓系数矩阵定 K：把 EERIE 的 4.635 分复算一遍

第三问的难点不在 K-means 本身，在于 K 到底取几。论文用了肘部法则和轮廓系数双重确认，这个做法值得抄，但很多人只画了肘部图就下结论——肘部图在 4 到 6 之间经常是平的，肉眼根本分不出来。

我的做法是直接算一个轮廓系数矩阵，把 K 从 2 到 8 全跑一遍，同时看惯性和轮廓系数两个指标：

from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler from sklearn.metrics import silhouette_score import pandas as pd Xs = StandardScaler().fit_transform(X) rows = [] for k in range(2, 9): km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(Xs) rows.append({ "k": k, "inertia": round(km.inertia_, 1), "silhouette": round(silhouette_score(Xs, km.labels_), 3), }) print(pd.DataFrame(rows))

n_init=10是必须的，K-means 对初始质心敏感，单次初始化可能落进局部最优；random_state=0保证结果可复现，写论文时别人照抄你的代码能跑出同样的数。

惯性看拐点，轮廓系数看峰值。两者冲突时以轮廓系数为准，因为它是「簇内紧、簇间远」的直接度量。这张表里 k=5 的轮廓系数最高，对应论文里分的五档难度，逻辑自洽。

定完 K，难度分的算法是：先用 PLS 预测该单词在七个尝试次数上的百分比，再按尝试次数加权求和。次数越高权重越大，七次以上（X）权重最高。EERIE 这个词复算出来是 4.635，落在第五档。它的特征组合也确实难：三个元音挨在一起、E重复三次、位置熵高，几个高难度信号叠在一起。

weights = np.array([1, 2, 3, 4, 5, 6, 8]) # 尝试次数权重，X 给最高 def difficulty_score(word): x = np.array([[build_features(word, freq, pos_dist, rank)[c] for c in X.columns]]) pct = pls.predict(x)[0] # 七个百分比 return float(np.dot(pct / pct.sum(), weights)) print(round(difficulty_score("eerie"), 3)) # 4.635 附近

权重这一列是可以调的。X给 8 而不是 7，是因为「七次以上仍未成功」在信息量上比「第七次成功」更强，给高一点能拉开档次。如果你的数据里 X 的占比普遍很低，把它调到 7 也行，但要保证五档之间的分界值重新算一遍。

最后一步验证，论文用的是混淆矩阵热图加准确率 0.941。混淆矩阵的对角线越集中，说明聚类结果和人工难度判断越一致。0.941 这个数在五分类里算相当高，但要注意前面 5.5 说的口径问题——这个数得是在独立验证集上算的，才站得住。

从那以后我每次做聚类题，都会先把轮廓系数矩阵打出来，再动手画肘部图，两个指标对不上就回去查标准化和特征口径，绝不凭一张拐点模糊的图拍板。希望帮到你。

本文还有配套的精品资源，点击获取