多因变量无标签小样本：ARIMA、PLS与K-means建模

简介：面向美赛参赛者与数据分析学习者的2023年美赛C题M论文，聚焦Wordle游戏数据的分析与预测，可作为竞赛论文写作、时间序列与多变量建模的学习参考，也适合毕业设计中的案例借鉴。全文围绕数据清洗与标准化、基于9个单词属性的Spearman相关性分析、时间序列区间预测、偏最小二乘回归以及K-means聚类展开，给出MAE约3～4、MSE约22～24的误差表现，并通过肘部法则与轮廓系数确定聚类数，结合难度评分和混淆矩阵评估模型，准确率达0.941。压缩包仅含1个PDF文件，大小约2.18MB，便于直接阅读与检索；目前已有170人学习。读者可从中获取完整赛题论文框架、问题拆解路径、模型选择依据、结果验证方式及不确定性讨论，对复盘美赛C题和迁移到类似数据分析任务具有较高参考价值。

1. 一道猜词游戏的数据题，为什么同时用上了 ARIMA、PLS 和 K-means

2023 年美赛 C 题给了一份 359 天的 Wordle 运营数据，最终要求回答四个问题：报告数量怎么变、某个词的七类尝试比例是多少、单词难度怎么分级、还有哪些有趣特征。表面看是数据分析预测，真正麻烦的地方在于：要预测的是一个七维百分比向量，每一维都受其他维约束，单词本身又没有任何现成难度标签。于是整套解法变成三层串接——ARIMA 管总量区间，PLS 管多因变量结构，K-means 管无标签分级。搞懂这三层为什么必须这么接，比记住任何单个模型都值钱，因为这类“多因变量 + 无标签 + 小样本”的组合，在毕业设计和工作里的用户行为数据上会反复出现。

2. 数据清洗与归一化：为什么 98% 和 102% 必须拉回 100%

2.1 两处断崖和两个四字母单词：清洗的第一现场

原始表有 359 行有效数据、12 列标签，缺失值为零，看起来干净得不像真实数据。但一查单词长度就露馅了：有两个词只有四个字母，对照游戏记录后修正为 clean 和 trash。这种错误不是笔误，而是数据从页面抓取或人工录入时被截断的典型痕迹，单词长度本身又是后续属性提取的输入，不修就会污染整条特征链。

真正影响时间序列的是两处断崖。2022/11/30 的报告数只有 2569，前后两天分别是 23739 和 22628；2022/2/13 的 Hard Mode 人数是 3249，前后是 9310 和 10343。处理方式是用两侧均值替代：报告数取 (23739+22628)/2≈23184，Hard Mode 取 (9310+10343)/2≈9827。这里没有选择直接删除这两天，原因很实在——ARIMA 要求等间隔时间序列，删掉一行会让时间索引错位，后面差分和滞后项全部对不上。

import pandas as pd import numpy as np df = pd.read_csv("wordle.csv") df["date"] = pd.to_datetime(df["date"]) df = df.sort_values("date").reset_index(drop=True) # 定位单点异常：与前后两日均值偏离超过 5 倍 for col in ["reported_results", "hard_mode_count"]: prev, nxt = df[col].shift(1), df[col].shift(-1) neighbor_mean = (prev + nxt) / 2 outlier = (df[col] < neighbor_mean * 0.2) | (df[col] > neighbor_mean * 5) # 只替换孤立异常点，连续异常不碰 df.loc[outlier, col] = neighbor_mean[outlier].round() # 七类百分比归一化到总和 100 pct_cols = ["try1", "try2", "try3", "try4", "try5", "try6", "tryX"] row_sum = df[pct_cols].sum(axis=1) df[pct_cols] = df[pct_cols].div(row_sum, axis=0) * 100

代码里的关键参数是偏离倍数阈值。取 5 倍是经验值：低于 2 倍容易把正常波动当异常，高于 10 倍又会漏掉 2569 这种只剩十分之一的塌陷。另一个细节是只替换孤立点，连续多天异常通常意味着统计口径变了，均值替代会掩盖真实结构。

2.2 均值替代的适用边界：什么时候不能这么干

均值替代在这份数据上成立，前提是异常点孤立、前后趋势平稳、且异常不携带信息。如果异常来自游戏改版、节日活动或统计口径调整，前后均值反而会把真实的结构变化抹平，让模型学到一个不存在的平滑曲线。常见的实现思路是先画时间序列图，看异常点是“尖刺”还是“台阶”：尖刺可以插值，台阶必须引入虚拟变量或分段建模。

这份数据的两处异常都是单日尖刺，前后没有持续偏移，所以均值替代是合理的。但即便合理，也应该在论文里记录替换前后的值，否则复现时没人知道 23184 这个数是怎么来的。血泪经验是：清洗步骤不写日志，后面调参时根本分不清误差是模型不行还是数据被改坏了。

2.3 百分比归一化的数学动机：单纯形约束

七类百分比相加，原始数据里出现了 98%、99%、100%、101%、102% 五种总和。来源不止一个：百分比四舍五入、分母用报告数还是完成数、X 类（七次及以上未成功）的统计边界。如果不处理，PLS 会把“行和不为 100”当成一个可学习的信号，回归系数就会去拟合这个人为误差，而不是单词属性与尝试分布的真实关系。

归一化的做法是每个百分比除以行和再乘 100，等价于把每个样本投影到七维空间里 Σp=100 的单纯形面上。这个约束很关键：七个因变量不是自由的，一个涨另一个必跌，协方差矩阵天然奇异，普通多元线性回归在这里会直接翻车。归一化之后，因变量的几何结构才和模型假设对齐。

提示：凡是遇到“若干比例加起来应该等于 100%”的数据，先跑一遍行和分布，再决定归一化还是补一个“其他”类，不要默认数据是干净的。

2.4 先可视化再建模：黄块分布说了什么

清洗完之后先画图，能得到几个建模前必须知道的先验。七类百分比堆叠图里，第四次尝试（黄色）占比最大，多数人集中在 3 到 5 次成功；未成功比例在 359 天里多数低于 5%，只有几十天高于 10%，说明出题方对难度的控制相当稳；Hard Mode 通关占比随时间逐步上升，说明玩家熟练度在积累。

这些结论看起来是描述统计，实际决定了后续模型的选择。比如未成功率长期低位，意味着 X 类百分比分布高度偏斜，直接对它做线性回归的残差不会好看；Hard Mode 占比有趋势，说明做时间序列时不能假设均值恒定。先可视化再建模，本质是给自己留一份“模型结果离谱时用来对照”的基准，否则黑匣子吐出什么数都只能照单全收。

3. ARIMA 预测：2023 年 3 月 1 日的报告数量区间是怎么算出来的

3.1 差分在消除什么：ARMA 与 ARIMA 的分界

ARMA(p,q) 的前提是序列平稳，均值、方差不随时间漂移。报告数量明显先升后降，直接喂 ARMA 会得到虚假回归。ARIMA 在 ARMA 前面加了一个 d 阶差分，把非平稳序列变成平稳序列：一阶差分是 y_t - y_{t-1}，二阶差分再对一阶差分做一次，写成算子形式就是 (1-B)^d y_t，B 是滞后算子。

差分不是越多越好。每差一次就丢一个样本，二阶差分还会放大噪声，让方差膨胀。判断 d 取多少，标准做法是单位根检验：ADF 检验在差分前不拒绝单位根，一阶差分后拒绝，就取 d=1；一阶差分后仍不拒绝再取 d=2。凭感觉定 d 是这类题里最常见的翻车点，差分过头会让预测区间虚胖，看起来“稳健”其实全是噪声。

3.2 p、d、q 定阶：ACF、PACF 与信息准则的配合

d 定下来之后，p 和 q 靠自相关函数和偏自相关函数的截尾、拖尾形态给初值，再用 AIC/BIC 在网格里挑。ACF 在 q 阶后截尾、PACF 拖尾，提示 MA(q)；PACF 在 p 阶后截尾、ACF 拖尾，提示 AR(p)；两者都拖尾就上 ARMA。实际流程通常是先网格搜索 p、q 从 0 到 5，按 AIC 排序，再对候选模型的残差做白噪声检验，残差还有自相关就继续加阶。

定阶阶段最容易被忽略的是残差检验。AIC 最小的模型不一定残差干净，如果残差里还留着周期性，说明模型没吃干净信息，预测区间会偏窄。常见的实现思路是画残差 ACF 图加 Ljung-Box 检验，p 值大于 0.05 才算过关。这一步不做，后面算出来的区间只是数学上成立，统计上站不住。

3.3 置信区间与预测区间：两个区间别混用

题目要的是 2023/3/1 报告数量的区间，必须是预测区间，不是置信区间。置信区间描述的是“均值响应”的不确定性，只包含参数估计误差；预测区间还要加上新观测本身的随机误差，所以更宽。公式上，预测区间标准差约等于 σ√(1 + x'(X'X)⁻¹x)，那个额外的 1 就是新观测的方差项。

很多人把回归输出的 confidence interval 直接当预测区间写进论文，区间偏窄，结论显得过于自信。区分方法很简单：问一句“这个区间是要覆盖一个未来的真实观测，还是覆盖未来观测的均值”，前者用预测区间，后者用置信区间。这道题显然是前者。

3.4 [17363,18994] 的边界条件

最终给出的 3 月 1 日预测区间是 [17363, 18994]，中值约 18178，区间宽度约 1600，相对宽度 9% 左右。这个区间成立有前提：历史模式延续、没有游戏改版、没有外部事件把流量打散、异常值处理没有改变趋势结构。3 月 1 日处于数据末端外推一步的位置，外推越远区间越宽，这是 ARIMA 的固有性质，不是模型缺陷。

读这个区间时要关注两件事。第一，区间是否覆盖了历史波动范围，如果历史同日波动远超这个宽度，说明模型低估了不确定性；第二，区间是否对异常值敏感，把 23184 换成原始 2569 再跑一遍，如果区间大幅偏移，说明清洗步骤对结论影响过大，需要在论文里做敏感性说明。这两条检查比区间本身更有信息量。

4. 从 PLS 到 K-means：七维因变量回归与无标签难度分级

4.1 七维因变量：为什么不能拆成七个独立回归

第二个问题要为每个词预测 (1,2,3,4,5,6,X) 七个百分比。直觉做法是跑七个线性回归，每个因变量一个模型。这个做法在这里有三重问题：七个百分比和为 100，因变量之间强负相关，七个模型的系数会互相矛盾；自变量是 9 个单词属性，属性之间本身共线（比如词频和首字母常见度）；样本只有 359 天，拆成七个模型后每个模型的有效信息更少。

PLS 的思路是把自变量和因变量同时投影到低维潜变量空间，在投影后的空间里做回归。它一次性处理全部七个因变量，潜变量提取时最大化 X 和 Y 的协方差，等于让模型自动找到“最能同时解释七个百分比”的方向。这比先 PCA 降维再回归多了一步监督信息，比直接多元回归多了一步降维，正好卡在这道题的需求上。

4.2 PLS 的潜变量提取与 MAE/MSE 的读法

PLS 的每次迭代做三件事：从 X 里提取一个与 Y 协方差最大的方向 w，算得分向量 t=Xw；对 Y 做同样的事得到 u；用 t 和 u 做回归并把残差传下去，重复提取下一个成分。成分数不能无限加，加多了会把噪声也拟合进去，通常用交叉验证的预测残差平方和定，取曲线拐点前的成分数。

论文给出的 MAE 在 3 到 4 之间，MSE 在 22 到 24 之间。注意单位是百分点，不是比例：平均每个百分比类别的预测误差约 3 到 4 个百分点，RMSE 约 4.7 到 4.9 个百分点。单类百分比常在 10% 到 30% 之间，5 个百分点的误差对应 20% 左右的相对误差，属于中等水平，不能算“高精度”。把 MAE 3.5 直接说成“误差很小”是过度解读，读这两个指标一定要带单位。

from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import cross_val_predict from sklearn.metrics import mean_absolute_error, mean_squared_error X = df[word_attr_cols].values # 9 个单词属性 Y = df[pct_cols].values # 7 类百分比，行和已归一化为 100 best = None for n in range(1, 8): # 成分数上限不超过自变量个数 pls = PLSRegression(n_components=n, scale=True) pred = cross_val_predict(pls, X, Y, cv=5) mae = mean_absolute_error(Y, pred) if best is None or mae < best[1]: best = (n, mae) print("最优成分数:", best[0], "MAE:", round(best[1], 2))

代码里 scale=True 是必须的，9 个属性量纲不同，不标准化会让数值大的属性主导协方差。成分数上限设为 7 而不是 9，是因为交叉验证在成分数接近自变量个数时容易过拟合，留一点余量更稳。这段代码跑完应该能复现 MAE 3 到 4 的量级，如果差很远，先检查因变量有没有归一化。

4.3 Spearman 相关：9 个属性的筛选与共线性

在进 PLS 之前，论文先对 9 个单词属性和 Hard Mode 得分百分比做了 Spearman 相关分析。选 Spearman 而不是 Pearson，是因为属性里有词频这类长尾分布，Pearson 会被极端值拉扯，Spearman 只看秩次，对单调非线性更稳。结论是 9 个属性对 Hard Mode 百分比都有一定影响，但“有一定影响”不等于“可以全部直接塞进线性回归”。

属性之间的共线性是真实存在的：首字母常见度和词频相关，字母重复度和单词长度相关。PLS 对共线性有一定容忍度，但相关分析的作用是给出解释框架——哪些属性是主要驱动，哪些是伴随效应。写结论时如果把 Spearman 显著直接说成因果，就会掉进“相关即因果”的坑，属性对难度的影响需要通过模型系数或变量投影重要性来交叉验证。

4.4 K-means 定 K：肘部法则的伪拐点与轮廓系数补位

第三个问题要给单词按难度分级，但数据里没有难度标签，只能无监督。K-means 的输入是单词属性加上 PLS 预测出的七类百分比，先标准化再聚类，否则量纲大的属性会垄断欧氏距离。K 的选择用了两个指标：肘部法则看簇内平方和随 K 下降的拐点，轮廓系数看每个点与自己簇的紧密度和与最近邻簇的分离度。

肘部法则的问题是拐点经常不唯一，K=3 和 K=5 都像拐点的情况很常见。轮廓系数 s=(b-a)/max(a,b)，a 是点到同簇其他点的平均距离，b 是到最近异簇的平均距离，取值在 -1 到 1 之间，越接近 1 分离度越好。两个指标一起看，再结合“五档难度”的业务可解释性，最终定 K=5。

from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score from sklearn.preprocessing import StandardScaler Xc = StandardScaler().fit_transform(cluster_features) for k in range(2, 9): km = KMeans(n_clusters=k, n_init=20, random_state=42).fit(Xc) s = silhouette_score(Xc, km.labels_) print(f"K={k}, SSE={km.inertia_:.1f}, Silhouette={s:.3f}")

n_init=20 是为了避开局部最优，K-means 对初始质心敏感，只跑一次容易得到不稳定的簇。SSE 用来画肘部曲线，Silhouette 用来交叉验证，两个都指向 5 附近才定 5。如果 SSE 拐点在 4、轮廓系数峰值在 6，说明数据结构本身模糊，这时候硬定 5 档就需要在论文里说明这是业务约束而不是数据结论。

4.5 难度加权评分与混淆矩阵：EERIE 的 4.635 和 0.941 准确率

定完簇之后，难度分不是簇标签直接给的，而是用加权公式算：把 PLS 预测的七类百分比按尝试次数赋权，尝试次数越多权重越高，X 类权重最高，加权求和得到连续难度分。EERIE 的得分是 4.635，落在第五档，说明这个词的预测分布偏向高尝试次数。这个分数比簇标签更有解释力，因为它保留了“差一点就进上一档”的信息。

模型评估用了混淆矩阵热图和准确率 0.941。这里要留个心眼：K-means 本身没有真实标签，混淆矩阵的“真值”来自哪里？如果是用同一批属性构造的代理标签，准确率就有循环论证的成分，0.941 衡量的是聚类结果与代理标签的一致性，不是与真实难度的吻合度。更稳的验证是留出一批人工标注的词做外部检验，或者在论文里明确说明代理标签的构造方式。0.941 可以写，但不能写成“模型难度分级准确率 94.1%”。

5. 避坑记录：五个让模型看起来很好看的陷阱

5.1 坑一：百分比不归一化就喂进模型

现象是 PLS 的 MAE 看起来很低，但预测出的七类百分比加起来是 103%，单类甚至出现负值。原因是行和不为 100 时，模型把总量误差也当成可学习的信号。解决办法是先做行和检查，把每个样本除以行和再乘 100，然后在预测端加一道非负裁剪和再归一化。这件事必须在建模前做，事后补救会让整个误差评估失真。

5.2 坑二：异常值一删了之

现象是删掉两天异常后 ARIMA 拟合得很好，但预测区间窄得离谱。原因是时间序列删点会破坏等间隔结构，滞后项错位，模型低估了真实波动。解决办法是孤立点用前后均值或插值替代，连续异常单独建虚拟变量，并在论文里记录替换前后的值。359 天的样本经不起删，每删一天都是信息损失。

5.3 坑三：只靠肘部法则定 K

现象是肘部曲线在 K=3 和 K=5 都有拐点，换一个随机种子聚类结果就变。原因是 SSE 单调下降，拐点判断带主观性。解决办法是肘部法则加轮廓系数交叉验证，再看每个簇的样本量和业务可解释性。如果三个指标打架，说明数据结构本身不适合硬分五档，应该在论文里承认这一点，而不是硬凑一个好看的 K。

5.4 坑四：0.941 的准确率直接写进结论

现象是混淆矩阵热图很漂亮，但每一类的召回率差距很大，小簇几乎全被吞进大簇。原因是准确率在类别不均衡时会被大簇主导。解决办法是看每一类的 precision、recall 和混淆矩阵对角线，而不是只看总准确率。如果第五档只有十几个词，0.941 里可能一个都没分对，这个数字就没有意义。

5.5 坑五：ARIMA 差分阶数凭感觉

现象是 d=2 之后残差看起来像白噪声，但预测区间比 d=1 宽了一倍。原因是对本来只需一阶差分的序列做了二次差分，过差分放大方差。解决办法是先跑 ADF 检验确定 d，差分后再看 ACF 是否还留趋势，残差做 Ljung-Box 检验。差分阶数是 ARIMA 里最容易被拍脑袋决定的参数，也是预测区间虚胖的主要来源。

6. 设计取舍：这套流程真正的价值不在 0.941 的准确率

三层模型的分工可以用一张表说清楚：

这张表暴露的核心取舍是：可解释性和预测精度在这里不是一回事。PLS 的 MAE 3 到 4 不算惊艳，但它把九个属性对七类百分比的影响拆开了，能回答“哪个属性让词变难”；K-means 的 0.941 看起来高，但它衡量的是聚类一致性，不是真实难度准确率。如果只追求预测精度，完全可以上梯度提升树，但那会丢掉属性层面的解释，而这道题第四问恰恰要的是“哪些属性与难度相关”。

加权公式里的权重是另一个主观点。尝试次数越多权重越高，这个方向没错，但具体给 1 次 1 分、X 类 7 分还是指数加权，会直接改变 EERIE 的 4.635 落在哪一档。更稳的做法是做权重敏感性分析：换两组权重，看排名前 10 的难词是否稳定。如果 EERIE 在两套权重下都进前五档，结论才站得住；如果只在某一套权重下成立，那 4.635 就只是一个数字，不是结论。

这套流程值得复用的地方，是它面对“多因变量 + 无标签 + 小样本”时的拆解顺序：先清洗到几何结构正确，再用 PLS 处理因变量之间的约束，最后用无监督方法补上缺失的标签层。这个顺序反过来做——先聚类再回归——会因为聚类输入里混进未归一化的百分比而整体跑偏。从那以后我每次拿到带百分比列的数据，都会先跑一遍行和检查，再决定要不要归一化；每次用 K-means 定 K，都会强制把轮廓系数和混淆矩阵一起打出来，不让单一指标替我做决定。希望帮到你。

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