Kaggle-工作手册-全-
Kaggle 工作手册（全）
原文：
annas-archive.org/md5/4c14de26a26e75bcd7b636b0f8e720d6译者：飞龙
第一章：01 最著名的表格竞赛：Porto Seguro 的 Safe Driver Prediction
加入我们的 Discord 书社区
packt.link/EarlyAccessCommunity
学习如何在任何 Kaggle 竞赛的排行榜上名列前茅需要耐心、勤奋和多次尝试，以便学会最佳的竞赛方式并取得优异成绩。因此，我们想到了一个工作簿，可以通过引导你尝试一些过去的 Kaggle 竞赛，并通过阅读讨论、重用笔记本、特征工程和训练各种模型来帮助你更快地掌握这些技能。
我们在书中从最著名的表格竞赛之一，Porto Seguro 的 Safe Driver Prediction 开始。在这个竞赛中，你需要解决保险中的一个常见问题，即确定下一年谁会提出汽车保险索赔。这样的信息对于提高那些更有可能提出索赔的司机的保险费，以及降低那些不太可能提出索赔的司机的保险费是有用的。
在阐述破解这个竞赛所需的关键洞察和技术细节时，我们将向您展示必要的代码，并要求您研究并回答在 Kaggle 书籍本身中可以找到的主题。因此，无需多言，让我们立即开始您的新学习之旅。
在本章中，你将学习：
-
如何调整和训练 LightGBM 模型。
-
如何构建去噪自编码器以及如何使用它来为神经网络提供数据。
-
如何有效地融合彼此差异很大的模型。
理解竞赛和数据
Porto Seguro 是巴西第三大保险公司（它在巴西和乌拉圭运营），提供汽车保险以及其他许多保险产品。在过去 20 年中，他们使用分析方法和机器学习来调整他们的价格，使汽车保险覆盖面更易于更多司机获得。为了探索实现他们任务的新方法，他们赞助了一个竞赛(www.kaggle.com/competitions/porto-seguro-safe-driver-prediction)，期望 Kagglers 能够提出解决他们核心分析问题的新方法。
该竞赛的目标是让 Kagglers 构建一个模型，预测司机在接下来一年内提出汽车保险索赔的概率，这是一个相当常见的任务（赞助商将其称为“保险的古典挑战”）。为此，赞助商提供了训练集和测试集，并且由于数据集不是很大，看起来准备得非常好，因此对任何人来说都是一个理想的竞赛。
如同在竞赛数据展示页面所述(www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/data)：
“属于相似分组特征在特征名称中标记为如此（例如，ind，reg，car，calc）。此外，特征名称包括后缀 bin 来指示二元特征，cat 来指示分类特征。没有这些指定的特征要么是连续的，要么是序数的。-1 的值表示该特征在观测中缺失。目标列表示该保单持有人是否提交了索赔”。
竞赛的数据准备非常仔细，以确保不泄露任何信息，尽管如此，关于特征的含义仍然保持保密，但很明显，可以将不同使用的标签指代到保险建模中常用的特定类型的特征：
-
ind 指代“个人特征”
-
“car”指的是“汽车特征”。
-
calc 指代“计算特征”
-
reg 指代“区域/地理特征”
至于个别特征，在比赛中已经有很多猜测。例如，参见www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/41489或www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/41488，这两篇文章都由 Raddar 撰写，或者再次尝试将特征归因于 Porto Seguro 的在线报价表www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/41057。尽管做出了所有这些努力，但直到现在，特征的含义仍然是个谜。
这个竞赛的有趣事实是：
-
数据是真实的，尽管特征是匿名的
-
数据准备得非常好，没有任何泄露（这里没有魔法特征）
-
测试集不仅与训练测试集具有相同的分类级别，而且似乎来自相同的分布，尽管山本优雅认为使用 t-SNE 对数据进行预处理会导致对抗验证测试失败（
www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/44784）。
作为第一次练习，参考 Kaggle Book 中关于对抗验证的内容和代码（从第 179 页开始），证明训练数据和测试数据很可能来自同一数据分布。
Tilii（蒙大拿州立大学的副教授 Mensur Dlakic）的一篇有趣的文章（www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/42197）展示了使用 tSNE 技术，指出“在保险参数方面，许多人非常相似，但其中一些人会提出索赔，而另一些人则不会”。Tilii 提到的情况在保险行业中很典型，即对于某些先验（保险参数）来说，发生某事的概率是相同的，但该事件是否发生则取决于我们观察情况的时间长短。
以保险行业中的物联网和遥测数据为例。分析驾驶员的行为以预测他们未来是否会提出索赔是很常见的。如果你的观察期太短（例如，像这次比赛一样，只有一年），那么即使是非常糟糕的驾驶员也可能不会提出索赔，因为这只是一个不太可能成为现实的低概率事件，而这一事件只有在一定时间后才会发生。类似的观点在 Andy Harless 的讨论中也有提及（www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/42735），他反而认为比赛的真正任务是猜测“一个潜在连续变量的值，该变量决定了哪些驾驶员更有可能发生事故”，因为实际上“提出索赔并不是驾驶员的特征；它是偶然的结果”。
理解评估指标
比赛中使用的指标是“标准化基尼系数”（以经济学中使用的类似基尼系数/指数命名），它之前在另一个比赛中也被使用过，即全美保险公司索赔预测挑战赛（www.kaggle.com/competitions/ClaimPredictionChallenge）。从那次比赛，我们可以清楚地了解这个指标的含义：
当你提交一个条目时，观察结果是从“预测值最大”到“预测值最小”排序的。这是唯一一个你的预测发挥作用的地方，所以只有由你的预测决定的顺序才是重要的。将观察结果从左到右可视化，预测值最大的在左边。然后我们从左到右移动，问“在数据的最左侧 x%中，你积累了多少实际观察到的损失？”如果没有模型，你可以预期在 10%的预测中积累 10%的损失，所以没有任何模型（或“零”模型）会得到一条直线。我们将你的曲线和这条直线之间的区域称为基尼系数。
“一个‘完美’模型有一个可达到的最大面积。我们将通过将你的模型的基尼系数除以完美模型的基尼系数来使用归一化基尼系数。”
竞赛中 Kilian Batzner 的笔记本提供了另一个很好的解释：www.kaggle.com/code/batzner/gini-coefficient-an-intuitive-explanation，尽管通过图表和玩具示例试图给一个不太常见的度量一个感觉，但在保险公司的精算部门。
在 Kaggle 书籍的第五章（第 95 页及以后），我们解释了如何处理比赛度量，特别是如果它们是新的且普遍未知的情况。作为一个练习，你能找出在 Kaggle 上有多少比赛使用了归一化基尼系数作为评估指标吗？
该度量可以通过 Mann–Whitney U 非参数统计检验和 ROC-AUC 分数来近似，因为它大约对应于 2 * ROC-AUC - 1。因此，最大化 ROC-AUC 等同于最大化归一化基尼系数（参见维基百科条目中的“与其他统计度量之间的关系”：en.wikipedia.org/wiki/Gini_coefficient）。
该度量也可以近似表示为缩放预测排名与缩放目标值的协方差，从而得到一个更易于理解的排名关联度量（参见 Dmitriy Guller：www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/40576）
从目标函数的角度来看，你可以优化二进制对数损失（就像你在分类问题中做的那样）。ROC-AUC 和归一化基尼系数都不是可微的，它们只能用于验证集上的度量评估（例如，用于早期停止或降低神经网络中的学习率）。然而，优化对数损失并不总是能提高 ROC-AUC 和归一化基尼系数。实际上，存在一个可微的 ROC-AUC 近似（Calders, Toon, 和 Szymon Jaroszewicz. "Efficient AUC optimization for classification." European conference on principles of data mining and knowledge discovery. Springer, Berlin, Heidelberg, 2007 link.springer.com/content/pdf/10.1007/978-3-540-74976-9_8.pdf）。然而，似乎在比赛中没有必要使用与对数损失不同的目标函数，以及将 ROC-AUC 或归一化基尼系数作为评估指标。
在这些笔记本中实际上有几个 Python 实现。我们在这里使用了 CPMP 的工作（www.kaggle.com/code/cpmpml/extremely-fast-gini-computation/notebook），该工作使用 Numba 来加速计算：它既精确又快速。
检查 Michael Jahrer 的顶级解决方案的想法
Michael Jahrer (www.kaggle.com/mjahrer，竞赛大师级选手，同时也是“BellKor's Pragmatic Chaos”团队在 Netflix Prize 竞赛中的获奖者之一)，在竞赛期间长时间以较大优势领先公开排行榜，并在最终私人排行榜公布后，被宣布为获胜者。
稍后，在讨论论坛中，他发布了他解决方案的简要总结，由于他对去噪自编码器和神经网络的巧妙使用，这个总结已成为许多 Kagglers 的参考（www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/44629）。尽管 Michael 没有附上关于他解决方案的任何 Python 代码（他将其称为“老式”和“低级”的，是直接用 C++/CUDA 编写的，没有使用 Python），但他的写作中充满了对所使用模型的引用，以及它们的超参数和架构。
首先，Michael 解释说，他的解决方案由六个模型（一个 LightGBM 模型和五个神经网络）的组合构成。此外，由于加权每个模型对组合的贡献（以及进行线性和非线性堆叠）可能无法获得优势，这可能是由于过拟合，他声称他转而使用了一个简单的模型组合（算术平均值），这些模型是从不同的种子构建的。
这样的洞察使我们的任务更容易复制他的方法，也因为他提到，仅仅将 LightGBM 的结果与他自己构建的神经网络的其中一个结果混合在一起，就足以保证在竞赛中获得第一名。这将限制我们的练习工作只关注两个优秀的单一模型，而不是一大堆模型。此外，他还提到，他做了很少的数据处理，但删除了一些列并对分类特征进行了独热编码。
构建 LightGBM 提交
我们的练习从基于 LightGBM 制定解决方案开始。您可以在以下地址找到已设置好的用于执行的代码：www.kaggle.com/code/lucamassaron/workbook-lgb。尽管我们已经提供了代码，但我们建议您直接从书中键入或复制代码并逐个执行代码块，理解每行代码的作用，并且您还可以进一步个性化解决方案，使其表现更佳。
我们首先导入关键包（Numpy、Pandas、Optuna 用于超参数优化、LightGBM 和一些实用函数）。我们还定义了一个配置类并实例化了它。随着我们对代码的探索，我们将讨论配置类中定义的参数。在此需要强调的是，通过使用包含所有参数的类，您将更容易在代码中一致地修改它们。在比赛的紧张时刻，很容易忘记更新在代码中多处引用的参数，并且当参数分散在单元格和函数中时，设置参数总是很困难。配置类可以节省您大量的精力并避免在过程中犯错误。
import numpy as np
import pandas as pd
import optuna
import lightgbm as lgb
from path import Path
from sklearn.model_selection import StratifiedKFold
class Config:
input_path = Path('../input/porto-seguro-safe-driver-prediction')
optuna_lgb = False
n_estimators = 1500
early_stopping_round = 150
cv_folds = 5
random_state = 0
params = {'objective': 'binary',
'boosting_type': 'gbdt',
'learning_rate': 0.01,
'max_bin': 25,
'num_leaves': 31,
'min_child_samples': 1500,
'colsample_bytree': 0.7,
'subsample_freq': 1,
'subsample': 0.7,
'reg_alpha': 1.0,
'reg_lambda': 1.0,
'verbosity': 0,
'random_state': 0}
config = Config()
下一步需要导入训练集、测试集和样本提交数据集。通过使用 pandas csv 读取函数，我们还设置了上传数据框的索引为每个数据示例的标识符（即‘id’列）。
由于属于相似分组特征被标记（使用 ind、reg、car、calc 标签在其标签中）以及二元和分类特征易于定位（它们分别使用 bin 和 cat 标签在其标签中），我们可以枚举它们并将它们记录到列表中。
train = pd.read_csv(config.input_path / 'train.csv', index_col='id')
test = pd.read_csv(config.input_path / 'test.csv', index_col='id')
submission = pd.read_csv(config.input_path / 'sample_submission.csv', index_col='id')
calc_features = [feat for feat in train.columns if "_calc" in feat]
cat_features = [feat for feat in train.columns if "_cat" in feat]
然后，我们只提取目标（一个由 0 和 1 组成的二元目标）并将其从训练数据集中删除。
target = train["target"]
train = train.drop("target", axis="columns")
到目前为止，正如 Michael Jahrer 所指出的，我们可以删除 calc 特征。这个想法在比赛中反复出现 (www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/41970)，尤其是在笔记本中，一方面是因为经验上可以验证删除它们可以提高公共排行榜的分数，另一方面是因为它们在梯度提升模型中的表现不佳（它们的重要性总是低于平均水平）。我们可以争论，由于它们是工程特征，它们在其原始特征方面不包含新信息，但它们只是向包含它们的任何训练模型添加噪声。
train = train.drop(calc_features, axis="columns")
test = test.drop(calc_features, axis="columns")
在比赛中，Tilii 使用了 Boruta 进行特征消除测试 (
github.com/scikit-learn-contrib/boruta_py)。您可以在以下链接找到他的内核：www.kaggle.com/code/tilii7/boruta-feature-elimination/notebook。如您所查，Boruta 没有将 calc_feature 作为确认的特征。练习：根据 Kaggle 书籍第 220 页提供的建议（“使用特征重要性评估你的工作”），作为一个练习，为这次比赛编写自己的特征选择笔记本，并检查应该保留哪些特征以及应该丢弃哪些特征。
分类特征则采用独热编码。由于我们想要重新训练它们的标签，并且由于相同的级别在训练和测试集中都存在（这是 Porto Seguro 团队在两个数据集之间进行仔细的训练/测试分割的结果），我们不是使用常规的 Scikit-Learn OneHotEncoder（scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OneHotEncoder.html），而是使用 pandas 的 get_dummies 函数（pandas.pydata.org/docs/reference/api/pandas.get_dummies.html）。由于 pandas 函数可能会根据特征及其级别在训练集和测试集中不同而产生不同的编码，我们确保独热编码的结果对于两者都是相同的。
train = pd.get_dummies(train, columns=cat_features)
test = pd.get_dummies(test, columns=cat_features)
assert((train.columns==test.columns).all())
在对分类特征进行独热编码后，我们已经完成了所有数据处理。我们继续定义我们的评估指标，即正常化基尼系数，正如之前讨论的那样。由于我们打算使用 LightGBM 模型，我们必须添加一个合适的包装器（gini_lgb），以便将训练集和验证集的评估以 LightGBM 算法可以处理的形式返回（见：lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.Booster.html?highlight=higher_better#lightgbm.Booster.eval - “每个评估函数应接受两个参数：preds, eval_data，并返回(eval_name, eval_result, is_higher_better)或此类元组的列表”）。
from numba import jit
@jit
def eval_gini(y_true, y_pred):
y_true = np.asarray(y_true)
y_true = y_true[np.argsort(y_pred)]
ntrue = 0
gini = 0
delta = 0
n = len(y_true)
for i in range(n-1, -1, -1):
y_i = y_true[i]
ntrue += y_i
gini += y_i * delta
delta += 1 - y_i
gini = 1 - 2 * gini / (ntrue * (n - ntrue))
return gini
def gini_lgb(y_true, y_pred):
eval_name = 'normalized_gini_coef'
eval_result = eval_gini(y_true, y_pred)
is_higher_better = True
return eval_name, eval_result, is_higher_better
关于训练参数，我们发现迈克尔·雅赫尔在其帖子中建议的参数（www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/44629）效果极佳。如果您在配置类中将optuna_lgb标志设置为 True，您也可以尝试通过 optuna（[optuna.org/](https://optuna.org/））进行搜索来找到相同的参数或类似性能的参数。在这里，优化尝试根据训练数据的 5 折交叉验证测试来找到关键参数（如学习率和正则化参数）的最佳值。为了加快速度，我们考虑了在验证过程中的早期停止（我们知道，这实际上可能有利于更好地拟合验证折的一些值 - 一个好的替代方案可能是移除早期停止回调并保持固定的训练轮数）。
if config.optuna_lgb:
def objective(trial):
params = {
'learning_rate': trial.suggest_float("learning_rate", 0.01, 1.0),
'num_leaves': trial.suggest_int("num_leaves", 3, 255),
'min_child_samples': trial.suggest_int("min_child_samples",
3, 3000),
'colsample_bytree': trial.suggest_float("colsample_bytree",
0.1, 1.0),
'subsample_freq': trial.suggest_int("subsample_freq", 0, 10),
'subsample': trial.suggest_float("subsample", 0.1, 1.0),
'reg_alpha': trial.suggest_loguniform("reg_alpha", 1e-9, 10.0),
'reg_lambda': trial.suggest_loguniform("reg_lambda", 1e-9, 10.0),
}
score = list()
skf = StratifiedKFold(n_splits=config.cv_folds, shuffle=True,
random_state=config.random_state)
for train_idx, valid_idx in skf.split(train, target):
X_train = train.iloc[train_idx]
y_train = target.iloc[train_idx]
X_valid = train.iloc[valid_idx]
y_valid = target.iloc[valid_idx]
model = lgb.LGBMClassifier(**params,
n_estimators=1500,
early_stopping_round=150,
force_row_wise=True)
callbacks=[lgb.early_stopping(stopping_rounds=150,
verbose=False)]
model.fit(X_train, y_train,
eval_set=[(X_valid, y_valid)],
eval_metric=gini_lgb, callbacks=callbacks)
score.append(
model.best_score_['valid_0']['normalized_gini_coef'])
return np.mean(score)
study = optuna.create_study(direction="maximize")
study.optimize(objective, n_trials=300)
print("Best Gini Normalized Score", study.best_value)
print("Best parameters", study.best_params)
params = {'objective': 'binary',
'boosting_type': 'gbdt',
'verbosity': 0,
'random_state': 0}
params.update(study.best_params)
else:
params = config.params
在比赛中，Tilii 测试了使用 Boruta（github.com/scikit-learn-contrib/boruta_py）进行特征消除。你可以在他的核函数中找到：www.kaggle.com/code/tilii7/boruta-feature-elimination/notebook。正如你可以检查的那样，没有 calc_feature 被 Boruta 视为已确认的特征。
在 Kaggle 书中，我们解释了超参数优化（从第 241 页开始）并提供了 LightGBM 模型的一些关键超参数。作为一个练习，尝试通过减少或增加探索的参数以及尝试 Scikit-Learn 中的随机搜索或减半搜索等替代方法来改进超参数搜索（第 245-246 页）。
一旦我们有了最佳参数（或者我们简单地尝试 Jahrer 的参数），我们就可以开始训练和预测。根据最佳解决方案的建议，我们的策略是在每个交叉验证折上训练一个模型，并使用该折来对测试预测的平均值做出贡献。代码片段将生成测试预测和训练集上的出卷预测，这后来对于确定如何集成结果非常有用。
preds = np.zeros(len(test))
oof = np.zeros(len(train))
metric_evaluations = list()
skf = StratifiedKFold(n_splits=config.cv_folds, shuffle=True, random_state=config.random_state)
for idx, (train_idx, valid_idx) in enumerate(skf.split(train,
target)):
print(f"CV fold {idx}")
X_train, y_train = train.iloc[train_idx], target.iloc[train_idx]
X_valid, y_valid = train.iloc[valid_idx], target.iloc[valid_idx]
model = lgb.LGBMClassifier(**params,
n_estimators=config.n_estimators,
early_stopping_round=config.early_stopping_round,
force_row_wise=True)
callbacks=[lgb.early_stopping(stopping_rounds=150),
lgb.log_evaluation(period=100, show_stdv=False)]
model.fit(X_train, y_train,
eval_set=[(X_valid, y_valid)],
eval_metric=gini_lgb, callbacks=callbacks)
metric_evaluations.append(
model.best_score_['valid_0']['normalized_gini_coef'])
preds += (model.predict_proba(test,
num_iteration=model.best_iteration_)[:,1]
/ skf.n_splits)
oof[valid_idx] = model.predict_proba(X_valid,
num_iteration=model.best_iteration_)[:,1]
模型训练不应该花费太多时间。最终，你可以在交叉验证过程中报告获得的归一化基尼系数。
print(f"LightGBM CV normalized Gini coefficient:
{np.mean(metric_evaluations):0.3f}
({np.std(metric_evaluations):0.3f})")
结果相当鼓舞人心，因为平均分数是 0.289，而值的方差相当小。
LightGBM CV Gini Normalized Score: 0.289 (0.015)
剩下的工作是将出卷和测试预测保存为提交，并在公共和私人排行榜上验证结果。
submission['target'] = preds
submission.to_csv('lgb_submission.csv')
oofs = pd.DataFrame({'id':train_index, 'target':oof})
oofs.to_csv('dnn_oof.csv', index=False)
获得的公共分数应该在 0.28442 左右。相关的私人分数约为 0.29121，将你在最终排行榜上的位置排在第 30 位。这是一个相当好的结果，但我们仍然需要将其与不同的模型，一个神经网络，进行混合。
尽管 Michael Jahrer 在他的帖子中提到，对训练集进行 Bagging（即对训练数据进行多次自助抽样并基于这些自助样本来训练多个模型）应该会增加性能，但增加的幅度并不大。
设置去噪自动编码器和深度神经网络
下一步不是设置一个去噪自动编码器（DAE）和一个可以从中学习和预测的神经网络。你可以在以下笔记本中找到运行代码：www.kaggle.com/code/lucamassaron/workbook-dae。笔记本可以在 GPU 模式下运行（更快），但也可以通过一些轻微的修改在 CPU 上运行。
你可以在 Kaggle 书中了解更多关于去噪自动编码器在 Kaggle 比赛中应用的信息，在第 226 页及以后。
实际上，在比赛中使用 DAEs 重现 Michael Jahrer 方法的例子并不多，但我们从 OsciiArt 在另一个比赛中使用的工作 TensorFlow 实现中汲取了经验www.kaggle.com/code/osciiart/denoising-autoencoder。
在这里，我们首先导入所有必要的包，特别是 TensorFlow 和 Keras。由于我们将创建多个神经网络，我们指出 TensorFlow 不要使用所有可用的 GPU 内存，通过使用实验性的set_memory_growth命令。这将避免在过程中出现内存溢出问题。我们还记录了 Leaky Relu 激活作为自定义激活，这样我们就可以在 Keras 层中通过字符串来提及它。
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from path import Path
import gc
import optuna
from sklearn.model_selection import StratifiedKFold
from scipy.special import erfinv
import tensorflow as tf
gpus = tf.config.experimental.list_physical_devices('GPU')
for gpu in gpus:
tf.config.experimental.set_memory_growth(gpu, True)
from tensorflow import keras
from tensorflow.keras import backend as K
from tensorflow.keras.layers import Input, Dense, BatchNormalization, Dropout
from tensorflow.keras.models import Model, load_model
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
from tensorflow.keras.regularizers import l2
from tensorflow.keras.metrics import AUC
from tensorflow.keras.utils import get_custom_objects
from tensorflow.keras.layers import Activation, LeakyReLU
get_custom_objects().update({'leaky-relu': Activation(LeakyReLU(alpha=0.2))})
与我们创建多个神经网络而不耗尽内存的意图相关，我们还定义了一个简单的函数来清理 GPU 中的内存并移除不再需要的模型。
def gpu_cleanup(objects):
if objects:
del(objects)
K.clear_session()
gc.collect()
我们还重新配置了 Config 类，以便考虑与去噪自动编码器和神经网络相关的多个参数。正如之前关于 LightGBM 所述，将所有参数放在一个地方可以简化以一致方式修改它们。
class Config:
input_path = Path('../input/porto-seguro-safe-driver-prediction')
dae_batch_size = 128
dae_num_epoch = 50
dae_architecture = [1500, 1500, 1500]
reuse_autoencoder = False
batch_size = 128
num_epoch = 150
units = [64, 32]
input_dropout=0.06
dropout=0.08
regL2=0.09
activation='selu'
cv_folds = 5
nas = False
random_state = 0
config = Config()
如前所述，我们加载数据集并继续处理特征，通过移除计算特征和将分类特征进行 one-hot 编码。我们保留缺失的案例值为-1，正如 Michael Jahrer 在他的解决方案中指出的。
train = pd.read_csv(config.input_path / 'train.csv', index_col='id')
test = pd.read_csv(config.input_path / 'test.csv', index_col='id')
submission = pd.read_csv(config.input_path / 'sample_submission.csv', index_col='id')
calc_features = [feat for feat in train.columns if "_calc" in feat]
cat_features = [feat for feat in train.columns if "_cat" in feat]
target = train["target"]
train = train.drop("target", axis="columns")
train = train.drop(calc_features, axis="columns")
test = test.drop(calc_features, axis="columns")
train = pd.get_dummies(train, columns=cat_features)
test = pd.get_dummies(test, columns=cat_features)
assert((train.columns==test.columns).all())
然而，与我们的先前方法不同，我们必须重新缩放所有非二进制或非 one-hot 编码的分类特征。重新缩放将允许自动编码器和神经网络优化算法更快地收敛到良好的解决方案，因为它将不得不在一个可比较和预定义的范围内处理值。与使用统计归一化不同，GaussRank 是一种允许将转换变量的分布修改为高斯分布的程序。
如某些论文所述，例如在批量归一化论文arxiv.org/pdf/1502.03167.pdf中，如果提供高斯输入，神经网络的表现会更好。根据这篇 NVIDIA 博客文章developer.nvidia.com/blog/gauss-rank-transformation-is-100x-faster-with-rapids-and-cupy/，GaussRank 在大多数情况下都有效，但当特征已经呈正态分布或极端不对称时（在这种情况下应用转换可能会降低性能）。
print("Applying GaussRank to columns: ", end='')
to_normalize = list()
for k, col in enumerate(train.columns):
if '_bin' not in col and '_cat' not in col and '_missing' not in col:
to_normalize.append(col)
print(to_normalize)
def to_gauss(x): return np.sqrt(2) * erfinv(x)
def normalize(data, norm_cols):
n = data.shape[0]
for col in norm_cols:
sorted_idx = data[col].sort_values().index.tolist()
uniform = np.linspace(start=-0.99, stop=0.99, num=n)
normal = to_gauss(uniform)
normalized_col = pd.Series(index=sorted_idx, data=normal)
data[col] = normalized_col
return data
train = normalize(train, to_normalize)
test = normalize(test, to_normalize)
我们可以在数据集的所有数值特征上分别应用 GaussRank 转换到训练和测试特征：
Applying GaussRank to columns: ['ps_ind_01', 'ps_ind_03', 'ps_ind_14', 'ps_ind_15', 'ps_reg_01', 'ps_reg_02', 'ps_reg_03', 'ps_car_11', 'ps_car_12', 'ps_car_13', 'ps_car_14', 'ps_car_15']
在归一化特征时，我们只需将我们的数据转换为 float32 类型的 NumPy 数组，这是 GPU 的理想输入。
features = train.columns
train_index = train.index
test_index = test.index
train = train.values.astype(np.float32)
test = test.values.astype(np.float32)
接下来，我们只准备一些有用的函数，如评估函数、归一化基尼系数和有助于表示 Keras 模型在训练集和验证集上拟合历史的绘图函数。
def plot_keras_history(history, measures):
rows = len(measures) // 2 + len(measures) % 2
fig, panels = plt.subplots(rows, 2, figsize=(15, 5))
plt.subplots_adjust(top = 0.99, bottom=0.01,
hspace=0.4, wspace=0.2)
try:
panels = [item for sublist in panels for item in sublist]
except:
pass
for k, measure in enumerate(measures):
panel = panels[k]
panel.set_title(measure + ' history')
panel.plot(history.epoch, history.history[measure],
label="Train "+measure)
try:
panel.plot(history.epoch,
history.history["val_"+measure],
label="Validation "+measure)
except:
pass
panel.set(xlabel='epochs', ylabel=measure)
panel.legend()
plt.show(fig)
from numba import jit
@jit
def eval_gini(y_true, y_pred):
y_true = np.asarray(y_true)
y_true = y_true[np.argsort(y_pred)]
ntrue = 0
gini = 0
delta = 0
n = len(y_true)
for i in range(n-1, -1, -1):
y_i = y_true[i]
ntrue += y_i
gini += y_i * delta
delta += 1 - y_i
gini = 1 - 2 * gini / (ntrue * (n - ntrue))
return gini
接下来的函数实际上要复杂一些，并且与去噪自动编码器和监督神经网络的功能更相关。batch_generator是一个函数，它将创建一个生成器，提供基于批次大小的数据块。它实际上不是一个独立的生成器，而是作为我们将要描述的更复杂的批次生成器的一部分，即mixup_generator。
def batch_generator(x, batch_size, shuffle=True, random_state=None):
batch_index = 0
n = x.shape[0]
while True:
if batch_index == 0:
index_array = np.arange(n)
if shuffle:
np.random.seed(seed=random_state)
index_array = np.random.permutation(n)
current_index = (batch_index * batch_size) % n
if n >= current_index + batch_size:
current_batch_size = batch_size
batch_index += 1
else:
current_batch_size = n - current_index
batch_index = 0
batch = x[index_array[current_index: current_index + current_batch_size]]
yield batch
mixup_generator是另一个生成器，它返回部分交换值的批次数据，以创建一些噪声并增强数据，以便去噪自动编码器（DAE）不会过度拟合训练集。它基于一个交换率，固定为 Michael Jahrer 建议的 15%的特征。
该函数生成两组不同的数据批次，一组用于释放给模型，另一组则用作在释放批次中交换值的来源。基于随机选择，其基本概率为交换率，在每个批次中，决定交换两个批次之间的一定数量的特征。
这样可以确保去噪自动编码器不能总是依赖于相同的特征（因为它们可能随时被随机交换），但它必须关注所有特征（在某种程度上类似于 dropout）。以便在它们之间找到关系，并在过程结束时正确地重建数据。
def mixup_generator(X, batch_size, swaprate=0.15, shuffle=True, random_state=None):
if random_state is None:
random_state = np.randint(0, 999)
num_features = X.shape[1]
num_swaps = int(num_features * swaprate)
generator_a = batch_generator(X, batch_size, shuffle,
random_state)
generator_b = batch_generator(X, batch_size, shuffle,
random_state + 1)
while True:
batch = next(generator_a)
mixed_batch = batch.copy()
effective_batch_size = batch.shape[0]
alternative_batch = next(generator_b)
assert((batch != alternative_batch).any())
for i in range(effective_batch_size):
swap_idx = np.random.choice(num_features, num_swaps,
replace=False)
mixed_batch[i, swap_idx] = alternative_batch[i, swap_idx]
yield (mixed_batch, batch)
get_DAE函数用于构建去噪自动编码器。它接受一个参数来定义架构，在我们的情况下，已经设置为三个各含 1500 个节点的层（如迈克尔·耶雷尔的建议）。第一层应作为编码器，第二层是瓶颈层，理想情况下包含能够表达数据信息的潜在特征，最后一层是解码层，能够重建初始输入数据。这三层具有 relu 激活函数，没有偏差，每一层后面都跟着一个批量归一化层。最终输出与重建的输入数据具有线性激活。训练使用具有标准设置的 adam 优化器（优化的成本函数是均方误差 - mse）。
def get_DAE(X, architecture=[1500, 1500, 1500]):
features = X.shape[1]
inputs = Input((features,))
for i, nodes in enumerate(architecture):
layer = Dense(nodes, activation='relu',
use_bias=False, name=f"code_{i+1}")
if i==0:
x = layer(inputs)
else:
x = layer(x)
x = BatchNormalization()(x)
outputs = Dense(features, activation='linear')(x)
model = Model(inputs=inputs, outputs=outputs)
model.compile(optimizer='adam', loss='mse',
metrics=['mse', 'mae'])
return model
这里仅报告extract_dae_features函数，仅用于教育目的。该函数有助于提取训练去噪自动编码器特定层的值。提取是通过构建一个新的模型来完成的，该模型结合了 DAE 输入层和所需的输出层。然后，简单的预测将提取我们需要的值（预测还允许固定首选的批次大小，以便满足任何内存需求）。
在竞赛的情况下，考虑到观察的数量和从自编码器中提取的特征数量，如果我们使用这个函数，得到的密集矩阵将太大，无法由 Kaggle 笔记本的内存处理。因此，我们的策略不会是将原始数据转换成瓶颈层的自编码器节点值，而是将自编码器与其冻结的层（直到瓶颈层）与监督神经网络融合，正如我们很快将要讨论的。
def extract_dae_features(autoencoder, X, layers=[3]):
data = []
for layer in layers:
if layer==0:
data.append(X)
else:
get_layer_output = Model([autoencoder.layers[0].input],
[autoencoder.layers[layer].output])
layer_output = get_layer_output.predict(X,
batch_size=128)
data.append(layer_output)
data = np.hstack(data)
return data
为了完成与 DAE 的工作，我们有一个最终函数，将所有之前的函数包装成一个无监督的训练过程（至少部分是无监督的，因为有一个早期停止监控器设置在验证集上）。该函数设置混合增强生成器，创建去噪自编码器架构，然后对其进行训练，监控其在验证集上的拟合度以实现早期停止，如果有过度拟合的迹象。最后，在返回训练好的 DAE 之前，它绘制了训练和验证拟合的图表，并将模型存储在磁盘上。
即使我们在该模型上尝试固定一个种子，与 LightGBM 模型相反，结果极其不稳定，它们可能会影响最终的集成结果。虽然结果可能是一个高分，但它可能会在私有排行榜上得分更高或更低，尽管在公共排行榜上获得的结果与公共排行榜非常相关，这将使你能够根据其公共结果始终选择最佳的最终提交。
def autoencoder_fitting(X_train, X_valid, filename='dae',
random_state=None, suppress_output=False):
if suppress_output:
verbose = 0
else:
verbose = 2
print("Fitting a denoising autoencoder")
tf.random.set_seed(seed=random_state)
generator = mixup_generator(X_train,
batch_size=config.dae_batch_size,
swaprate=0.15,
random_state=config.random_state)
dae = get_DAE(X_train, architecture=config.dae_architecture)
steps_per_epoch = np.ceil(X_train.shape[0] /
config.dae_batch_size)
early_stopping = EarlyStopping(monitor='val_mse',
mode='min',
patience=5,
restore_best_weights=True,
verbose=0)
history = dae.fit(generator,
steps_per_epoch=steps_per_epoch,
epochs=config.dae_num_epoch,
validation_data=(X_valid, X_valid),
callbacks=[early_stopping],
verbose=verbose)
if not suppress_output: plot_keras_history(history,
measures=['mse', 'mae'])
dae.save(filename)
return dae
在处理了 DAE 之后，我们也有机会定义一个监督神经网络模型，该模型应该预测我们的索赔预期。作为第一步，我们定义了一个函数来定义工作中的一个单层：
-
随机正态初始化，因为经验上已经发现这种方法在这个问题中能收敛到更好的结果
-
一个具有 L2 正则化和可参数化激活函数的密集层
-
可排除和可调整的 dropout
下面是创建密集块的代码：
def dense_blocks(x, units, activation, regL2, dropout):
kernel_initializer = keras.initializers.RandomNormal(mean=0.0,
stddev=0.1, seed=config.random_state)
for k, layer_units in enumerate(units):
if regL2 > 0:
x = Dense(layer_units, activation=activation,
kernel_initializer=kernel_initializer,
kernel_regularizer=l2(regL2))(x)
else:
x = Dense(layer_units,
kernel_initializer=kernel_initializer,
activation=activation)(x)
if dropout > 0:
x = Dropout(dropout)(x)
return x
如您可能已经注意到的，定义单层的函数相当可定制。同样，对于包装架构函数也是如此，它接受层的数量和其中的单元数、dropout 概率、正则化和激活类型作为输入。我们的想法是能够运行神经架构搜索（NAS），并找出在我们的问题中应该表现更好的配置。
关于该函数的最后一句话，在输入中需要提供训练好的 DAE，因为它的输入被用作神经网络模型的输入，而它的第一层与 DAE 的输出相连。这样，我们实际上是将两个模型合并为一个（DAE 的权重已经冻结，不可训练）。这个解决方案是为了避免需要转换所有训练数据，而只需要神经网络处理的单个批次，从而在系统中节省内存。
def dnn_model(dae, units=[4500, 1000, 1000],
input_dropout=0.1, dropout=0.5,
regL2=0.05,
activation='relu'):
inputs = dae.get_layer("code_2").output
if input_dropout > 0:
x = Dropout(input_dropout)(inputs)
else:
x = tf.keras.layers.Layer()(inputs)
x = dense_blocks(x, units, activation, regL2, dropout)
outputs = Dense(1, activation='sigmoid')(x)
model = Model(inputs=dae.input, outputs=outputs)
model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001),
loss=keras.losses.binary_crossentropy,
metrics=[AUC(name='auc')])
return model
我们以一个训练过程的包装来总结，包括训练整个管道所需的步骤，以交叉验证折进行训练。
def model_fitting(X_train, y_train, X_valid, y_valid, autoencoder,
filename, random_state=None, suppress_output=False):
if suppress_output:
verbose = 0
else:
verbose = 2
print("Fitting model")
early_stopping = EarlyStopping(monitor='val_auc',
mode='max',
patience=10,
restore_best_weights=True,
verbose=0)
rlrop = ReduceLROnPlateau(monitor='val_auc',
mode='max',
patience=2,
factor=0.75,
verbose=0)
tf.random.set_seed(seed=random_state)
model = dnn_model(autoencoder,
units=config.units,
input_dropout=config.input_dropout,
dropout=config.dropout,
regL2=config.regL2,
activation=config.activation)
history = model.fit(X_train, y_train,
epochs=config.num_epoch,
batch_size=config.batch_size,
validation_data=(X_valid, y_valid),
callbacks=[early_stopping, rlrop],
shuffle=True,
verbose=verbose)
model.save(filename)
if not suppress_output:
plot_keras_history(history, measures=['loss', 'auc'])
return model, history
由于我们的 DAE 实现与 Jahrer 的不同，尽管背后的想法相同，我们不能完全依赖他对监督神经网络架构的指示，我们必须寻找理想的架构，就像我们在 LightGBM 模型中寻找最佳超参数一样。使用 Optuna 并利用我们为配置网络架构而设置的多个参数，我们可以运行这个代码片段几个小时，并了解什么可能工作得更好。
在我们的实验中，我们注意到：
-
我们应该使用一个具有少量节点的两层网络，分别是 64 和 32。
-
输入 dropout、层间 dropout 以及一些 L2 正则化确实有所帮助。
-
使用 SELU 激活函数会更好。
这里是运行整个优化实验的代码片段：
if config.nas is True:
def evaluate():
metric_evaluations = list()
skf = StratifiedKFold(n_splits=config.cv_folds, shuffle=True, random_state=config.random_state)
for k, (train_idx, valid_idx) in enumerate(skf.split(train, target)):
X_train, y_train = train[train_idx, :], target[train_idx]
X_valid, y_valid = train[valid_idx, :], target[valid_idx]
if config.reuse_autoencoder:
autoencoder = load_model(f"./dae_fold_{k}")
else:
autoencoder = autoencoder_fitting(X_train, X_valid,
filename=f'./dae_fold_{k}',
random_state=config.random_state,
suppress_output=True)
model, _ = model_fitting(X_train, y_train, X_valid, y_valid,
autoencoder=autoencoder,
filename=f"dnn_model_fold_{k}",
random_state=config.random_state,
suppress_output=True)
val_preds = model.predict(X_valid, batch_size=128, verbose=0)
best_score = eval_gini(y_true=y_valid, y_pred=np.ravel(val_preds))
metric_evaluations.append(best_score)
gpu_cleanup([autoencoder, model])
return np.mean(metric_evaluations)
def objective(trial):
params = {
'first_layer': trial.suggest_categorical("first_layer", [8, 16, 32, 64, 128, 256, 512]),
'second_layer': trial.suggest_categorical("second_layer", [0, 8, 16, 32, 64, 128, 256]),
'third_layer': trial.suggest_categorical("third_layer", [0, 8, 16, 32, 64, 128, 256]),
'input_dropout': trial.suggest_float("input_dropout", 0.0, 0.5),
'dropout': trial.suggest_float("dropout", 0.0, 0.5),
'regL2': trial.suggest_uniform("regL2", 0.0, 0.1),
'activation': trial.suggest_categorical("activation", ['relu', 'leaky-relu', 'selu'])
}
config.units = [nodes for nodes in [params['first_layer'], params['second_layer'], params['third_layer']] if nodes > 0]
config.input_dropout = params['input_dropout']
config.dropout = params['dropout']
config.regL2 = params['regL2']
config.activation = params['activation']
return evaluate()
study = optuna.create_study(direction="maximize")
study.optimize(objective, n_trials=60)
print("Best Gini Normalized Score", study.best_value)
print("Best parameters", study.best_params)
config.units = [nodes for nodes in [study.best_params['first_layer'], study.best_params['second_layer'], study.best_params['third_layer']] if nodes > 0]
config.input_dropout = study.best_params['input_dropout']
config.dropout = study.best_params['dropout']
config.regL2 = study.best_params['regL2']
config.activation = study.best_params['activation']
如果你想要了解更多关于神经网络架构搜索（NAS）的信息，你可以查看 Kaggle 书，从第 276 页开始。在 DAE 和监督神经网络的情况下，寻找最佳架构至关重要，因为我们正在实施与 Michael Jahrer 解决方案肯定不同的东西。
作为一项练习，尝试通过使用 KerasTuner（可在 Kaggle 书中第 285 页及以后找到）来改进超参数搜索，这是一种优化神经网络的快速解决方案，它得到了 Keras 的创造者 François Chollet 的重要贡献。
在一切准备就绪之后，我们就可以开始训练了。大约一个小时后，在带有 GPU 的 Kaggle 笔记本上，你可以获得完整的测试和跨折预测。
preds = np.zeros(len(test))
oof = np.zeros(len(train))
metric_evaluations = list()
skf = StratifiedKFold(n_splits=config.cv_folds, shuffle=True, random_state=config.random_state)
for k, (train_idx, valid_idx) in enumerate(skf.split(train, target)):
print(f"CV fold {k}")
X_train, y_train = train[train_idx, :], target[train_idx]
X_valid, y_valid = train[valid_idx, :], target[valid_idx]
if config.reuse_autoencoder:
print("restoring previously trained dae")
autoencoder = load_model(f"./dae_fold_{k}")
else:
autoencoder = autoencoder_fitting(X_train, X_valid,
filename=f'./dae_fold_{k}',
random_state=config.random_state)
model, history = model_fitting(X_train, y_train, X_valid, y_valid,
autoencoder=autoencoder,
filename=f"dnn_model_fold_{k}",
random_state=config.random_state)
val_preds = model.predict(X_valid, batch_size=128)
best_score = eval_gini(y_true=y_valid,
y_pred=np.ravel(val_preds))
best_epoch = np.argmax(history.history['val_auc']) + 1
print(f"[best epoch is {best_epoch}]\tvalidation_0-gini_dnn: {best_score:0.5f}\n")
metric_evaluations.append(best_score)
preds += (model.predict(test, batch_size=128).ravel() /
skf.n_splits)
oof[valid_idx] = model.predict(X_valid, batch_size=128).ravel()
gpu_cleanup([autoencoder, model])
就像我们对 LighGBM 模型所做的那样，我们可以通过查看平均折归一化基尼系数来了解结果。
print(f"DNN CV normalized Gini coefficient: {np.mean(metric_evaluations):0.3f} ({np.std(metric_evaluations):0.3f})")
结果不会与之前使用 LightGBM 获得的结果完全一致。
DNN CV Gini Normalized Score: 0.276 (0.015)
制作提交并提交将导致公开分数大约为 0.27737，私人分数大约为 0.28471（结果可能与我们之前提到的有很大差异），并不是一个很高的分数。
submission['target'] = preds
submission.to_csv('dnn_submission.csv')
oofs = pd.DataFrame({'id':train_index, 'target':oof})
oofs.to_csv('dnn_oof.csv', index=False)
神经网络的少量结果似乎遵循了这样一个谚语：神经网络在表格问题上的表现不佳。无论如何，作为 Kagglers，我们知道所有模型都对在排行榜上取得成功有用，我们只需要找出如何最好地使用它们。当然，一个使用自动编码器的神经网络已经提出了一种受数据噪声影响较小、以不同方式阐述信息的解决方案，这比 GBM 要好。
结果集成
现在，拥有两个模型后，剩下的就是将它们混合在一起，看看我们是否能提高结果。正如 Jahrer 所建议的，我们直接尝试将它们混合，但我们并不局限于仅仅产生两个模型的平均值（因为我们的方法最终与 Jahrer 的方法略有不同），我们还将尝试为混合找到最佳权重。我们开始导入折叠外的预测，并准备好我们的评估函数。
import pandas as pd
import numpy as np
from numba import jit
@jit
def eval_gini(y_true, y_pred):
y_true = np.asarray(y_true)
y_true = y_true[np.argsort(y_pred)]
ntrue = 0
gini = 0
delta = 0
n = len(y_true)
for i in range(n-1, -1, -1):
y_i = y_true[i]
ntrue += y_i
gini += y_i * delta
delta += 1 - y_i
gini = 1 - 2 * gini / (ntrue * (n - ntrue))
return gini
lgb_oof = pd.read_csv("../input/workbook-lgb/lgb_oof.csv")
dnn_oof = pd.read_csv("../input/workbook-dae/dnn_oof.csv")
target = pd.read_csv("../input/porto-seguro-safe-driver-prediction/train.csv", usecols=['id','target'])
一旦完成，我们将 LightGBM 和神经网络的折叠外预测转换为排名，因为归一化的基尼系数对排名敏感（就像 ROC-AUC 评估一样）。
lgb_oof_ranks = (lgb_oof.target.rank() / len(lgb_oof))
dnn_oof_ranks = (dnn_oof.target.rank() / len(dnn_oof))
现在我们只是测试，通过使用不同的权重结合两个模型，我们是否能得到更好的折叠外数据评估。
baseline = eval_gini(y_true=target.target, y_pred=lgb_oof_ranks)
print(f"starting from a oof lgb baseline {baseline:0.5f}\n")
best_alpha = 1.0
for alpha in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]:
ensemble = alpha * lgb_oof_ranks + (1.0 - alpha) * dnn_oof_ranks
score = eval_gini(y_true=target.target, y_pred=ensemble)
print(f"lgd={alpha:0.1f} dnn={(1.0 - alpha):0.1f} -> {score:0.5f}")
if score > baseline:
baseline = score
best_alpha = alpha
print(f"\nBest alpha is {best_alpha:0.1f}")
准备就绪后，通过运行代码片段，我们可以得到有趣的结果：
starting from a oof lgb baseline 0.28850
lgd=0.1 dnn=0.9 -> 0.27352
lgd=0.2 dnn=0.8 -> 0.27744
lgd=0.3 dnn=0.7 -> 0.28084
lgd=0.4 dnn=0.6 -> 0.28368
lgd=0.5 dnn=0.5 -> 0.28595
lgd=0.6 dnn=0.4 -> 0.28763
lgd=0.7 dnn=0.3 -> 0.28873
lgd=0.8 dnn=0.2 -> 0.28923
lgd=0.9 dnn=0.1 -> 0.28916
Best alpha is 0.8
看起来，使用强权重（0.8）在 LightGBM 模型上和弱权重（0.2）在神经网络上混合可能会产生一个表现优异的模型。我们立即通过为模型设置相同的权重和我们已经找到的理想权重来测试这个假设。
lgb_submission = pd.read_csv("../input/workbook-lgb/lgb_submission.csv")
dnn_submission = pd.read_csv("../input/workbook-dae/dnn_submission.csv")
submission = pd.read_csv(
"../input/porto-seguro-safe-driver-prediction/sample_submission.csv")
首先，我们尝试等权重解决方案：
lgb_ranks = (lgb_submission.target.rank() / len(lgb_submission))
dnn_ranks = (dnn_submission.target.rank() / len(dnn_submission))
submission.target = lgb_ranks * 0.5 + dnn_ranks * 0.5
submission.to_csv("equal_blend_rank.csv", index=False)
这导致了公共得分为 0.28393，私人得分为 0.29093，大约位于最终排行榜的第 50 位，离我们的预期有点远。现在让我们尝试使用折叠外预测帮助我们找到的权重：
lgb_ranks = (lgb_submission.target.rank() / len(lgb_submission))
dnn_ranks = (dnn_submission.target.rank() / len(dnn_submission))
submission.target = lgb_ranks * best_alpha + dnn_ranks * (1.0 - best_alpha)
submission.to_csv("blend_rank.csv", index=False)
这里，结果导致公共得分为 0.28502，私人得分为 0.29192，最终在最终排行榜上大约位于第七位。确实是一个更好的结果，因为 LightGBM 是一个很好的模型，但它可能缺少一些数据中的细微差别，这些差别可以通过添加来自在去噪数据上训练的神经网络的某些信息来作为有利的纠正。
如 CPMP 在他的解决方案中指出的（
www.kaggle.com/competitions/porto-seguro-safe-driver-prediction/discussion/44614），根据如何构建你的交叉验证，你可能会经历“折叠间基尼分数的巨大变化”。因此，CPMP 建议通过使用多个交叉验证的许多不同种子来减少估计的方差，并平均结果。练习：作为一个练习，尝试修改我们使用的代码，以便创建更稳定的预测，特别是对于去噪自编码器。
摘要
在本章中，你已经处理了一个经典的表格竞赛。通过阅读竞赛的笔记本和讨论，我们提出了一种简单的解决方案，仅涉及两个易于混合的模型。特别是，我们提供了一个示例，说明如何使用去噪自编码器来生成一种特别适用于处理表格数据的替代数据处理方法。通过理解和复制过去竞赛中的解决方案，你可以在 Kaggle 竞赛中快速建立你的核心能力，并迅速在最近的竞赛和挑战中表现出色和稳定。
在下一章中，我们将探索 Kaggle 的另一个表格竞赛，这次是关于一个复杂的时间序列预测问题。
第二章：2 Makridakis 竞赛：M5 在 Kaggle 上的准确性和不确定性
加入我们的 Discord 书籍社区
packt.link/EarlyAccessCommunity
自 1982 年以来，Spyros Makridakis (mofc.unic.ac.cy/dr-spyros-makridakis/https://mofc.unic.ac.cy/dr-spyros-makridakis/) 一直让来自世界各地的研究团队参与预测挑战，称为 M 竞赛，以便对现有和新预测方法在不同预测问题上的有效性进行比较。因此，M 竞赛始终对学术界和实践者完全开放。这些竞赛可能是预测社区中被引用和参考最多的活动，并且它们始终突出了预测方法领域的最新进展。每个先前的 M 竞赛都为研究人员和实践者提供了有用的数据来训练和测试他们的预测工具，同时也提供了一系列发现和途径，这些发现和途径正在彻底改变预测的方式。
最近举行的 M5 竞赛（M6 竞赛正在撰写本章时进行）在 Kaggle 上举行，它证明了在尝试解决一系列零售产品销量预测问题时，梯度提升方法的有用性尤为重要。在本章中，我们专注于准确性赛道，处理了 Kaggle 竞赛中的一个时间序列问题，通过复制一个顶级、简单且最清晰的解决方案之一，我们旨在为读者提供代码和思路，以成功应对未来可能出现在 Kaggle 上的任何预测竞赛。
-
除了竞赛页面，我们在《国际预测杂志》的以下论文中找到了关于竞赛及其动态的大量信息：
-
Makridakis, Spyros, Evangelos Spiliotis, 和 Vassilios Assimakopoulos. 《M5 竞赛：背景、组织和实施》。《国际预测杂志》（2021 年）。
-
Makridakis, Spyros, Evangelos Spiliotis, 和 Vassilios Assimakopoulos。"M5 准确性竞赛：结果、发现和结论"《国际预测杂志》（2022 年）。
-
Makridakis, Spyros, 等人。"《M5 不确定性竞赛：结果、发现和结论》"《国际预测杂志》（2021 年）。
理解竞争和数据
竞赛从 2020 年 3 月持续到 6 月，超过 7,000 名参与者参加了在 Kaggle 上的竞赛。组织者将其分为两个独立的赛道，一个用于点预测（准确性赛道），另一个用于在不同置信区间估计可靠值（不确定性赛道）。
沃尔玛提供了数据。它包括 42,840 个每日销售时间序列，这些序列按部门、类别和店铺分层排列，分布在三个美国州（这些序列彼此之间有一定的相关性）。除了销售数据，沃尔玛还提供了伴随信息（外生变量，通常在预测问题中不常提供），例如商品价格、一些日历信息、相关的促销活动或其他影响销售的事件。
除了 Kaggle，数据以及之前 M 竞赛的数据集都可以在以下地址找到：forecasters.org/resources/time-series-data/。
竞赛的一个有趣方面是它处理了快速消费品和慢速消费品销售，有许多例子展示了最新的间歇性销售（销售通常是零，但在一些罕见情况下）。虽然间歇性序列在许多行业中很常见，但对于许多从业者来说，在预测中仍然是一个具有挑战性的案例。
竞赛时间表被安排为两部分。在第一部分，从 2020 年 3 月初到 6 月 1 日，参赛者可以在 1,913 天范围内的任何一天训练模型，并在公共测试集（从 1,914 天到 1,941 天）上对其提交进行评分。在那之后，直到 7 月 1 日的竞赛结束，公共测试集作为训练集的一部分提供，允许参与者调整他们的模型以预测从 1,942 天到 1969 天（28 天的时间窗口，即四周）。在那个时期，提交在排行榜上没有得分。
竞赛这种安排背后的比例是为了让团队最初能在排行榜上测试他们的模型，并为他们提供在笔记本和讨论中分享最佳表现方法的基础。在第一阶段之后，组织者希望避免排行榜被用于过度拟合或模型超参数调整的目的，他们希望模拟一种预测情况，就像在现实世界中发生的那样。此外，只选择一个提交作为最终提交的要求，反映了现实世界的相同必要性（在现实世界中，你不能使用两个不同的模型预测，然后在之后选择最适合你的一个）。
关于数据，我们提到数据由沃尔玛提供，并代表美国市场：它源自加利福尼亚州、威斯康星州和德克萨斯州的 10 家商店。具体来说，数据由 3,049 产品的销售组成，分为三个类别（爱好、食品和家庭），每个类别可以进一步分为 7 个部门。这种层次结构无疑是一个挑战，因为你可以在美国市场、州市场、单个商店、产品类别、类别部门和最终特定产品级别建模销售动态。所有这些级别也可以组合成不同的汇总，这是第二轨道，即不确定性轨道中需要预测的内容：
| 级别 ID | 级别描述 | 汇总级别 | 序列数量 |
|---|---|---|---|
| 1 | 所有产品，按所有商店和州汇总 | 总计 | 1 |
| 2 | 所有产品，按每个州汇总 | 州 | 3 |
| 3 | 所有产品，按每个商店汇总 | 商店 | 10 |
| 4 | 所有产品，按每个类别汇总 | 类别 | 3 |
| 5 | 所有产品，按每个部门汇总 | 部门 | 7 |
| 6 | 所有产品，按每个州和类别汇总 | 州-类别 | 9 |
| 7 | 所有产品，按每个州和部门汇总 | 州-部门 | 21 |
| 8 | 所有产品，按每个商店和类别汇总 | 商店-类别 | 30 |
| 9 | 所有产品，按每个商店和部门汇总 | 商店-部门 | 70 |
| 10 | 每个产品，按所有商店/州汇总 | 产品 | 3,049 |
| 11 | 每个产品，按每个州汇总 | 产品-州 | 9,147 |
| 12 | 每个产品，按每个商店汇总 | 产品-商店 | 30,490 |
| 总计 | 42,840 |
从时间角度来看，粒度是每日销售记录，覆盖了从 2011 年 1 月 29 日到 2016 年 6 月 19 日的期间，总计 1,969 天，其中 1,913 天用于训练，28 天用于验证 - 公开排行榜 - 28 天用于测试 - 私人排行榜。实际上，在零售行业中，28 天的预测范围被认为是处理大多数商品的库存和重新订购操作的正确范围。
让我们检查一下比赛中收到的不同数据。你将获得 sales_train_evaluation.csv、sell_prices.csv 和 calendar.csv。其中包含时间序列的是 sales_train_evaluation.csv。它由作为标识符的字段（item_id、dept_id、cat_id、store_id 和 state_id）以及从 d_1 到 d_1941 的列组成，代表那些天的销售情况：
图 2.1：sales_train_evaluation.csv 数据
sell_prices.csv 包含关于商品价格的信息。这里的难点在于将 wm_yr_wk（周标识符）与训练数据中的列连接起来：
图 2.2：sell_prices.csv 数据
最后一个文件，calendar.csv，包含可能影响销售的事件相关数据：
图 2.3：calendar.csv 数据
同样，主要困难似乎在于将数据与训练表中的列连接起来。无论如何，在这里您可以获得一个简单的键来连接列（d 字段）与wm_yr_wk。此外，在表中，我们表示了可能发生在特定日期的不同事件，以及 SNAP 日，这是特别的日子，在这些日子里，可以使用的营养援助福利补充营养援助计划（SNAP）。
理解评估指标
准确度竞赛引入了一个新的评估指标：加权均方根缩放误差（WRMSSE）。该指标评估点预测围绕预测序列实现值平均值的偏差：
其中：
n 是训练样本的长度
h 是预测范围（在我们的情况下是h=28）
Y[t] 是时间 t 的销售价值，
是时间 t 的预测值
在竞赛指南(mofc.unic.ac.cy/m5-competition/)中，关于 WRMSSE，指出：
-
RMSSE 的分母仅针对那些正在积极销售的考察产品的时间段进行计算，即，在评估序列观察到的第一个非零需求之后的时期。
-
该度量与规模无关，这意味着它可以有效地用于比较不同规模序列的预测。
-
与其他度量相比，它可以安全地计算，因为它不依赖于可能等于或接近零的值的除法（例如，在Y[t] = 0 时进行的百分比误差，或者当用于缩放的基准误差为零时进行的相对误差）。
-
该度量对正负预测误差、大预测和小预测进行同等惩罚，因此是对称的。
亚历山大·索阿雷（Alexander Soare）在这篇帖子中提供了对这种工作原理的良好解释(www.kaggle.com/alexandersoare)：www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/148273。在转换评估指标后，亚历山大将性能的提高归因于预测误差与日销售价值日变化率之间的比率提高。如果误差与日变化率相同（比率=1），则模型可能并不比基于历史变化的随机猜测好多少。如果您的误差小于这个值，则分数以二次方式（由于平方根）接近零。因此，WRMSSE 为 0.5 对应比率为 0.7，WRMSSE 为 0.25 对应比率为 0.5。
在竞赛期间，许多尝试不仅将指标用于排行榜的评估，还将其作为目标函数。首先，Tweedie 损失（在 XGBoost 和 LightGBM 中实现）对于这个问题来说效果相当好，因为它可以处理大多数产品的销售分布的偏斜（其中很多也有间歇性销售，这也被 Tweedie 损失很好地处理）。泊松和伽马分布可以被认为是 Tweedie 分布的极端情况：基于参数幂，p，当 p=1 时得到泊松分布，当 p=2 时得到伽马分布。这种幂参数实际上是连接分布的均值和方差的粘合剂，通过公式方差 = k*均值**p。使用介于 1 和 2 之间的幂值，实际上可以得到泊松和伽马分布的混合，这可以很好地适应竞赛问题。实际上，大多数参与竞赛并使用 GBM 解决方案的 Kagglers 都求助于 Tweedie 损失。
尽管 Tweedie 方法取得了成功，然而，一些其他 Kagglers 却发现了一些有趣的方法来实现一个更接近 WRMSSE 的目标损失，用于他们的模型：
-
Martin Kovacevic Buvinic 及其非对称损失：
www.kaggle.com/code/ragnar123/simple-lgbm-groupkfold-cv/notebook -
Timetraveller 使用 PyTorch Autograd 获取梯度和对角线，以实现任何可微分的连续损失函数，并在 LighGBM 中实现：
www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/152837
检验来自 Monsaraida 的第四名解决方案的想法
对于这次竞赛，有许多解决方案可供选择，主要可以在竞赛 Kaggle 讨论页面上找到。挑战的前五名方法也已被竞赛组织者自己收集并发布（但有一个因为版权问题）：github.com/Mcompetitions/M5-methods（顺便说一句，重现获奖提交的结果是收集竞赛奖金的前提条件）。
显然，所有在比赛中获得较高排名的 Kagglers 都使用了 LightGBM，作为他们独特的模型类型或在集成/堆叠中，因为其较低的内存使用量和计算速度，这使其在处理和预测大量时间序列时在竞争中具有优势。但还有其他原因导致其成功。与基于 ARIMA 的经典方法相反，它不需要依赖于自相关分析，以及在具体确定问题中每个单个序列的参数。此外，与基于深度学习的方法相反，它不需要寻找改进复杂的神经网络架构或调整大量超参数。梯度提升方法在时间序列问题中的优势（对于扩展其他梯度提升算法，例如 XGBoost）是依赖于特征工程，根据时间滞后、移动平均和序列属性分组创建正确的特征数量。然后选择正确的目标函数并进行一些超参数调整，当时间序列足够长时（对于较短的序列，ARIMA 或指数平滑等经典统计方法仍然是推荐的选择），就足以获得优秀的结果。
与比赛中深度学习解决方案相比，LightGBM 和 XGBoost 的另一个优势是 Tweedie 损失，不需要任何特征缩放（深度学习网络对所使用的缩放特别敏感）以及训练速度，这允许在测试特征工程时进行更快的迭代。
在所有这些可用的解决方案中，我们发现由日本计算机科学家 Monsaraida（Masanori Miyahara）提出的方案最为有趣。他提出了一种简单直接的方法，在私人排行榜上排名第 4，得分为 0.53583。该方案仅使用一般特征，没有进行预先选择（如销售统计、日历、价格和标识符）。此外，它使用有限数量的同类型模型，使用 LightGBM 梯度提升，没有求助于任何类型的混合、递归建模（当预测反馈给其他层次相关的预测或乘数时，即选择常数以更好地拟合测试集）。以下是他在 M（github.com/Mcompetitions/M5-methods/tree/master/Code%20of%20Winning%20Methods/A4）的演示解决方案中提出的方案，可以注意到他对待每个商店的每个四周，最终对应于产生 40 个模型：
图 2.4：Monsaraida 关于其解决方案结构的说明
由于 Monsaraida 保持了其解决方案的简单和实用，就像在现实世界的预测项目中一样，在这一章中，我们将尝试通过重构他的代码来复制他的示例，以便在 Kaggle 笔记本中运行（我们将通过将代码拆分成多个笔记本来处理内存和运行时间限制）。这样，我们旨在为读者提供一个简单而有效的方法，基于梯度提升，来处理预测问题。
计算特定日期和时间跨度的预测
复制 Monsaraida 解决方案的计划是创建一个可由输入参数自定义的笔记本，以便生成训练和测试所需的必要处理数据以及用于预测的 LightGBM 模型。给定过去的数据，这些模型将被训练以学习预测未来特定天数内的值。通过让每个模型学习预测未来特定周范围内的值，可以获得最佳结果。由于我们必须预测未来最多 28 天，我们需要一个从未来+1 天到未来+7 天的模型，然后是另一个能够从未来+8 天到未来+14 天的模型，再是另一个从未来+15 天到+21 天的模型，最后是一个能够处理从未来+22 天到未来+28 天的预测的模型。我们需要为每个时间范围创建一个 Kaggle 笔记本，因此我们需要四个笔记本。每个这些笔记本都将被训练以预测每个参与竞赛的十个店铺的未来时间跨度。总共，每个笔记本将生成十个模型。所有这些笔记本将共同生成四十个模型，覆盖所有未来范围和所有店铺。
由于我们需要为公共排行榜和私有排行榜都进行预测，因此有必要重复这个过程两次，在 1,913 天（预测 1,914 天到 1,941 天的日子）停止训练以提交公共测试集，以及在 1,941 天（预测 1,942 天到 1,969 天的日子）停止训练以提交私有测试集。
由于基于 CPU 运行 Kaggle 笔记本的限制，所有这八个笔记本都可以并行运行（整个过程几乎需要 6 个半小时）。每个笔记本可以通过其名称与其他笔记本区分开来，包含与最后训练日和前瞻性预测天数相关的参数值。这些笔记本中的一个示例可以在以下链接找到：www.kaggle.com/code/lucamassaron/m5-train-day-1941-horizon-7。
现在我们一起来检查代码是如何安排的，以及我们可以从 Monsaraida 的解决方案中学到什么。
我们首先导入必要的包。您只需注意到，除了 NumPy 和 pandas 之外，唯一的数据科学专用包是 LightGBM。您可能还会注意到，我们将使用 gc（垃圾回收）：这是因为我们需要限制脚本使用的内存量，并且我们经常只是收集和回收未使用的内存。作为这一策略的一部分，我们也经常将模型和数据结构存储到磁盘上，而不是保留在内存中：
import numpy as np
import pandas as pd
import os
import random
import math
from decimal import Decimal as dec
import datetime
import time
import gc
import lightgbm as lgb
import pickle
import warnings
warnings.filterwarnings(“ignore”, category=UserWarning)
作为限制内存使用的策略的一部分，我们求助于 Kaggle 书中描述的用于减少 pandas DataFrame 的函数，最初由 Arjan Groen 在 Zillow 比赛中开发（阅读讨论www.kaggle.com/competitions/tabular-playground-series-dec-2021/discussion/291844）：
def reduce_mem_usage(df, verbose=True):
numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']
start_mem = df.memory_usage().sum() / 1024**2
for col in df.columns:
col_type = df[col].dtypes
if col_type in numerics:
c_min = df[col].min()
c_max = df[col].max()
if str(col_type)[:3] == 'int':
if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
df[col] = df[col].astype(np.int8)
elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
df[col] = df[col].astype(np.int16)
elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
df[col] = df[col].astype(np.int32)
elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
df[col] = df[col].astype(np.int64)
else:
if c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
df[col] = df[col].astype(np.float32)
else:
df[col] = df[col].astype(np.float64)
end_mem = df.memory_usage().sum() / 1024**2
if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))
return df
我们继续定义这个解决方案的函数，因为这样有助于将解决方案拆分成更小的部分，并且当您从函数返回时，清理所有使用的变量会更加容易（您只需保留保存到磁盘的内容，并从函数返回）。我们的下一个函数帮助我们加载所有可用的数据并将其压缩：
def load_data():
train_df = reduce_mem_usage(pd.read_csv("../input/m5-forecasting-accuracy/sales_train_evaluation.csv"))
prices_df = reduce_mem_usage(pd.read_csv("../input/m5-forecasting-accuracy/sell_prices.csv"))
calendar_df = reduce_mem_usage(pd.read_csv("../input/m5-forecasting-accuracy/calendar.csv"))
submission_df = reduce_mem_usage(pd.read_csv("../input/m5-forecasting-accuracy/sample_submission.csv"))
return train_df, prices_df, calendar_df, submission_df
train_df, prices_df, calendar_df, submission_df = load_data()
在准备获取价格、体积和日历信息相关的数据代码之后，我们继续准备第一个处理函数，该函数将具有创建一个基本的信息表的角色，其中item_id、dept_id、cat_id、state_id和store_id作为行键，一个日期列和一个包含体积的值列。这是通过使用 pandas 命令 melt(pandas.pydata.org/pandas-docs/stable/reference/api/pandas.melt.html)从具有所有日期数据列的行开始的。该命令以 DataFrame 的索引为参考，然后选择所有剩余的特征，将它们的名称放在一个列上，将它们的值放在另一个列上（var_name和value_name参数帮助您定义这些新列的名称）。这样，您可以将代表某个商店中某个商品的销售额序列的行展开成多个行，每行代表一天。展开列的顺序保持不变，这保证了现在您的时间序列在垂直轴上（因此您可以在其上应用进一步的转换，如移动平均值）。
为了让您了解正在发生的情况，这里展示了在pd.melt转换之前的train_df。注意，不同日期的量被作为列特征：
图 2.5：训练数据框
变换之后，您将获得一个grid_df，其中日期被分配到单独的日期上：
图 2.5：将 pd.melt 应用于训练数据框
特征 d 包含了对不属于索引的列的引用，本质上，是从 d_1 到 d_1935 的所有特征。通过简单地从其值中移除‘d_’前缀并将它们转换为整数，你现在就有一个日特征。
除了这个之外，代码片段还根据时间从训练数据中分离出一部分行（你的验证集）。在训练部分，它还会根据你提供的预测范围添加必要的行来进行预测。
下面是创建我们的基本特征模板的函数。作为输入，它接受 train_df DataFrame，它期望训练结束的日期和预测范围（你想要预测的未来天数）：
def generate_base_grid(train_df, end_train_day_x, predict_horizon):
index_columns = ['id', 'item_id', 'dept_id', 'cat_id', 'store_id', 'state_id']
grid_df = pd.melt(train_df, id_vars=index_columns, var_name='d', value_name='sales')
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df['d_org'] = grid_df['d']
grid_df['d'] = grid_df['d'].apply(lambda x: x[2:]).astype(np.int16)
time_mask = (grid_df['d'] > end_train_day_x) & (grid_df['d'] <= end_train_day_x + predict_horizon)
holdout_df = grid_df.loc[time_mask, ["id", "d", "sales"]].reset_index(drop=True)
holdout_df.to_feather(f"holdout_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(holdout_df)
gc.collect()
grid_df = grid_df[grid_df['d'] <= end_train_day_x]
grid_df['d'] = grid_df['d_org']
grid_df = grid_df.drop('d_org', axis=1)
add_grid = pd.DataFrame()
for i in range(predict_horizon):
temp_df = train_df[index_columns]
temp_df = temp_df.drop_duplicates()
temp_df['d'] = 'd_' + str(end_train_day_x + i + 1)
temp_df['sales'] = np.nan
add_grid = pd.concat([add_grid, temp_df])
grid_df = pd.concat([grid_df, add_grid])
grid_df = grid_df.reset_index(drop=True)
for col in index_columns:
grid_df[col] = grid_df[col].astype('category')
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
在处理完创建基本特征模板的函数之后，我们准备了一个 pandas DataFrame 的合并函数，这有助于在处理大量数据时节省内存空间并避免内存错误。给定两个 DataFrame，df1 和 df2 以及我们需要它们合并的外键集合，该函数在 df1 和 df2 之间应用左外连接，而不创建新的合并对象，只是简单地扩展现有的 df1 DataFrame。
该函数首先从 df1 中提取外键，然后将提取的键与 df2 合并。这样，该函数创建了一个新的 DataFrame，称为 merged_gf，其顺序与 df1 相同。在这个时候，我们只是将 merged_gf 列分配给 df1。内部，df1 将从 merged_gf 中选择对内部数据结构的引用。这种做法有助于最小化内存使用，因为在任何时刻都只创建了必要使用的数据（没有可以填充内存的重复数据）。当函数返回 df1 时，merged_gf 被取消，但 df1 现在使用的数据。
下面是这个实用函数的代码：
def merge_by_concat(df1, df2, merge_on):
merged_gf = df1[merge_on]
merged_gf = merged_gf.merge(df2, on=merge_on, how='left')
new_columns = [col for col in list(merged_gf)
if col not in merge_on]
df1[new_columns] = merged_gf[new_columns]
return df1
在完成这个必要的步骤之后，我们继续编写一个新的数据处理函数。这次我们处理价格数据，这是一组包含每个店铺每个商品所有周的价格的数据集。由于确定我们是否在谈论一个新商品出现在店铺中非常重要，该函数选择了价格可用的第一个日期（使用价格表中的 wm_yr_wk 特征，代表周 id）并将其复制到我们的特征模板中。
下面是处理发布日期的代码：
def calc_release_week(prices_df, end_train_day_x, predict_horizon):
index_columns = ['id', 'item_id', 'dept_id', 'cat_id', 'store_id', 'state_id']
grid_df = pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
release_df = prices_df.groupby(['store_id', 'item_id'])['wm_yr_wk'].agg(['min']).reset_index()
release_df.columns = ['store_id', 'item_id', 'release']
grid_df = merge_by_concat(grid_df, release_df, ['store_id', 'item_id'])
del release_df
grid_df = reduce_mem_usage(grid_df, verbose=False)
gc.collect()
grid_df = merge_by_concat(grid_df, calendar_df[['wm_yr_wk', 'd']], ['d'])
grid_df = grid_df.reset_index(drop=True)
grid_df['release'] = grid_df['release'] - grid_df['release'].min()
grid_df['release'] = grid_df['release'].astype(np.int16)
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
在处理完产品在店铺中出现的日期之后，我们肯定要继续处理价格。就每个商品而言，每个店铺，我们准备了一些基本的价格特征，告诉我们：
-
实际价格（按最大值归一化）
-
最高价格
-
最低价格
-
平均价格
-
价格的标准差
-
该商品所经历的不同价格数量
-
店铺中具有相同价格的物品数量
除了这些基本的价格描述性统计之外，我们还添加了一些特征来描述每个商品在不同时间粒度下的动态变化：
-
日 momentum，即实际价格与其前一天价格的比例
-
月 momentum，即实际价格与其同月平均价格的比例
-
年 momentum，即实际价格与其同年平均价格的比例
在这里，我们使用了两个有趣且必要的 pandas 方法进行时间序列特征处理：
-
shift：可以将索引向前或向后移动 n 步（
pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.shift.html） -
transform：应用于 group by，用转换后的值填充类似索引的特征（
pandas.pydata.org/docs/reference/api/pandas.core.groupby.DataFrameGroupBy.transform.html）
此外，为了揭示商品以心理定价阈值（例如$19.99 或£2.98——参见此讨论：www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/145011）出售的情况，处理了价格的小数部分作为特征。math.modf函数（docs.python.org/3.8/library/math.html#math.modf）有助于这样做，因为它将任何浮点数分成小数部分和整数部分（一个两项元组）。
最后，结果表被保存到磁盘上。
下面是这个函数对价格进行所有特征工程的过程：
def generate_grid_price(prices_df, calendar_df, end_train_day_x, predict_horizon):
grid_df = pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
prices_df['price_max'] = prices_df.groupby(['store_id', 'item_id'])['sell_price'].transform('max')
prices_df['price_min'] = prices_df.groupby(['store_id', 'item_id'])['sell_price'].transform('min')
prices_df['price_std'] = prices_df.groupby(['store_id', 'item_id'])['sell_price'].transform('std')
prices_df['price_mean'] = prices_df.groupby(['store_id', 'item_id'])['sell_price'].transform('mean')
prices_df['price_norm'] = prices_df['sell_price'] / prices_df['price_max']
prices_df['price_nunique'] = prices_df.groupby(['store_id', 'item_id'])['sell_price'].transform('nunique')
prices_df['item_nunique'] = prices_df.groupby(['store_id', 'sell_price'])['item_id'].transform('nunique')
calendar_prices = calendar_df[['wm_yr_wk', 'month', 'year']]
calendar_prices = calendar_prices.drop_duplicates(subset=['wm_yr_wk'])
prices_df = prices_df.merge(calendar_prices[['wm_yr_wk', 'month', 'year']], on=['wm_yr_wk'], how='left')
del calendar_prices
gc.collect()
prices_df['price_momentum'] = prices_df['sell_price'] / prices_df.groupby(['store_id', 'item_id'])[
'sell_price'].transform(lambda x: x.shift(1))
prices_df['price_momentum_m'] = prices_df['sell_price'] / prices_df.groupby(['store_id', 'item_id', 'month'])[
'sell_price'].transform('mean')
prices_df['price_momentum_y'] = prices_df['sell_price'] / prices_df.groupby(['store_id', 'item_id', 'year'])[
'sell_price'].transform('mean')
prices_df['sell_price_cent'] = [math.modf(p)[0] for p in prices_df['sell_price']]
prices_df['price_max_cent'] = [math.modf(p)[0] for p in prices_df['price_max']]
prices_df['price_min_cent'] = [math.modf(p)[0] for p in prices_df['price_min']]
del prices_df['month'], prices_df['year']
prices_df = reduce_mem_usage(prices_df, verbose=False)
gc.collect()
original_columns = list(grid_df)
grid_df = grid_df.merge(prices_df, on=['store_id', 'item_id', 'wm_yr_wk'], how='left')
del(prices_df)
gc.collect()
keep_columns = [col for col in list(grid_df) if col not in original_columns]
grid_df = grid_df[['id', 'd'] + keep_columns]
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"grid_price_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
下一个函数计算月相，返回其八个阶段之一（从新月到亏月）。尽管月相不应直接影响任何销售（天气条件反而会，但我们没有数据中的天气信息），但它们代表了一个 29 天半的周期，这非常适合周期性购物行为。关于为什么月相可能作为预测因素的不同假设，在这次竞赛帖子中有有趣的讨论：www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/154776：
def get_moon_phase(d): # 0=new, 4=full; 4 days/phase
diff = datetime.datetime.strptime(d, '%Y-%m-%d') - datetime.datetime(2001, 1, 1)
days = dec(diff.days) + (dec(diff.seconds) / dec(86400))
lunations = dec("0.20439731") + (days * dec("0.03386319269"))
phase_index = math.floor((lunations % dec(1) * dec(8)) + dec('0.5'))
return int(phase_index) & 7
月相函数是创建基于时间特征的通用函数的一部分。该函数接受日历数据集信息并将其放置在特征中。此类信息包含事件及其类型，以及指示 SNAP 时期（一种名为补充营养援助计划的营养援助福利——简称 SNAP——以帮助低收入家庭）的信息，这些信息可能进一步推动基本商品的销售。该函数还生成诸如日期、月份、年份、星期几、月份中的星期以及是否为周末等数值特征。以下是代码：
def generate_grid_calendar(calendar_df, end_train_day_x, predict_horizon):
grid_df = pd.read_feather(
f"grid_df_{end_train_day_x}_to_{end_train_day_x +
predict_horizon}.feather")
grid_df = grid_df[['id', 'd']]
gc.collect()
calendar_df['moon'] = calendar_df.date.apply(get_moon_phase)
# Merge calendar partly
icols = ['date',
'd',
'event_name_1',
'event_type_1',
'event_name_2',
'event_type_2',
'snap_CA',
'snap_TX',
'snap_WI',
'moon',
]
grid_df = grid_df.merge(calendar_df[icols], on=['d'], how='left')
icols = ['event_name_1',
'event_type_1',
'event_name_2',
'event_type_2',
'snap_CA',
'snap_TX',
'snap_WI']
for col in icols:
grid_df[col] = grid_df[col].astype('category')
grid_df['date'] = pd.to_datetime(grid_df['date'])
grid_df['tm_d'] = grid_df['date'].dt.day.astype(np.int8)
grid_df['tm_w'] = grid_df['date'].dt.isocalendar().week.astype(np.int8)
grid_df['tm_m'] = grid_df['date'].dt.month.astype(np.int8)
grid_df['tm_y'] = grid_df['date'].dt.year
grid_df['tm_y'] = (grid_df['tm_y'] - grid_df['tm_y'].min()).astype(np.int8)
grid_df['tm_wm'] = grid_df['tm_d'].apply(lambda x: math.ceil(x / 7)).astype(np.int8)
grid_df['tm_dw'] = grid_df['date'].dt.dayofweek.astype(np.int8)
grid_df['tm_w_end'] = (grid_df['tm_dw'] >= 5).astype(np.int8)
del(grid_df['date'])
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"grid_calendar_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
del(calendar_df)
gc.collect()
以下函数只是移除了wm_yr_wk特征，并将 d（日）特征转换为数值。这是以下特征转换函数的必要步骤。
def modify_grid_base(end_train_day_x, predict_horizon):
grid_df = pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
grid_df['d'] = grid_df['d'].apply(lambda x: x[2:]).astype(np.int16)
del grid_df['wm_yr_wk']
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
我们最后两个特征创建函数将为时间序列生成更复杂的特征工程。第一个函数将生成滞后销售和它们的移动平均值。首先，使用移动方法（pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.shift.html）将生成过去 15 天的滞后销售范围。然后使用移动方法结合滚动（pandas.pydata.org/docs/reference/api/pandas.DataFrame.rolling.html）将创建 7 天、14 天、30 天、60 天和 180 天的移动平均值。
移动命令是必要的，因为它将允许移动索引，这样你就可以始终考虑可用于计算的数据。因此，如果你的预测范围达到七天，计算将只考虑七天前的数据。然后滚动命令将创建一个移动窗口观察值，可以总结（在这种情况下是通过平均值）。在一个时期（移动窗口）上有一个平均值，并跟踪其演变，将帮助你更好地检测趋势中的任何变化，因为不会在时间窗口中重复的图案将被平滑。这是时间序列分析中去除噪声和非有趣模式的一种常见策略。例如，使用七天的滚动平均值，你可以取消所有日间模式，仅表示你的销售在每周发生的情况。
你可以尝试不同的移动平均值窗口吗？尝试不同的策略可能会有所帮助。例如，通过探索 2022 年 1 月的 Tabular Playground 竞赛（
www.kaggle.com/competitions/tabular-playground-series-jan-2022），该竞赛致力于时间序列，你可能会找到更多的想法，因为大多数解决方案都是使用梯度提升构建的。
下面是生成滞后和滚动平均值特征的代码：
def generate_lag_feature(end_train_day_x, predict_horizon):
grid_df = pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
grid_df = grid_df[['id', 'd', 'sales']]
num_lag_day_list = []
num_lag_day = 15
for col in range(predict_horizon, predict_horizon + num_lag_day):
num_lag_day_list.append(col)
grid_df = grid_df.assign(**{
'{}_lag_{}'.format(col, l): grid_df.groupby(['id'])['sales'].transform(lambda x: x.shift(l))
for l in num_lag_day_list
})
for col in list(grid_df):
if 'lag' in col:
grid_df[col] = grid_df[col].astype(np.float16)
num_rolling_day_list = [7, 14, 30, 60, 180]
for num_rolling_day in num_rolling_day_list:
grid_df['rolling_mean_' + str(num_rolling_day)] = grid_df.groupby(['id'])['sales'].transform(
lambda x: x.shift(predict_horizon).rolling(num_rolling_day).mean()).astype(np.float16)
grid_df['rolling_std_' + str(num_rolling_day)] = grid_df.groupby(['id'])['sales'].transform(
lambda x: x.shift(predict_horizon).rolling(num_rolling_day).std()).astype(np.float16)
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"lag_feature_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
至于第二个高级特征工程函数，它是一个编码函数，它接受变量在州、商店、类别、部门和销售商品中的特定分组，并代表它们的平均值和标准差。这种嵌入是时间无关的（时间不是分组的一部分），并且它们的作用是帮助训练算法区分商品、类别和商店（及其组合）如何相互区分。
如同《Kaggle 书》第 216 页所述，使用目标编码计算所提出的嵌入相当简单，你能获得更好的结果吗？
代码通过分组特征，计算它们的描述性统计（在我们的情况下是均值或标准差），然后使用我们之前讨论过的 transform(pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.transform.html)方法将结果应用于数据集：
def generate_target_encoding_feature(end_train_day_x, predict_horizon):
grid_df = pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
grid_df.loc[grid_df['d'] > (end_train_day_x - predict_horizon), 'sales'] = np.nan
base_cols = list(grid_df)
icols = [
['state_id'],
['store_id'],
['cat_id'],
['dept_id'],
['state_id', 'cat_id'],
['state_id', 'dept_id'],
['store_id', 'cat_id'],
['store_id', 'dept_id'],
['item_id'],
['item_id', 'state_id'],
['item_id', 'store_id']
]
for col in icols:
col_name = '_' + '_'.join(col) + '_'
grid_df['enc' + col_name + 'mean'] = grid_df.groupby(col)['sales'].transform('mean').astype(
np.float16)
grid_df['enc' + col_name + 'std'] = grid_df.groupby(col)['sales'].transform('std').astype(
np.float16)
keep_cols = [col for col in list(grid_df) if col not in base_cols]
grid_df = grid_df[['id', 'd'] + keep_cols]
grid_df = reduce_mem_usage(grid_df, verbose=False)
grid_df.to_feather(f"target_encoding_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
完成了特征工程部分后，我们现在继续将存储在磁盘上的所有文件组合起来，同时生成特征。以下函数仅加载基本特征、价格特征、日历特征、滞后/滚动和嵌入特征的不同数据集，并将它们全部连接起来。然后代码仅过滤与特定商店相关的行，将其保存为单独的数据集。这种做法与针对特定商店训练模型以预测特定时间间隔的策略相匹配：
def assemble_grid_by_store(train_df, end_train_day_x, predict_horizon):
grid_df = pd.concat([pd.read_feather(f"grid_df_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather"),
pd.read_feather(f"grid_price_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather").iloc[:, 2:],
pd.read_feather(f"grid_calendar_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather").iloc[:, 2:]],
axis=1)
gc.collect()
store_id_set_list = list(train_df['store_id'].unique())
index_store = dict()
for store_id in store_id_set_list:
extract = grid_df[grid_df['store_id'] == store_id]
index_store[store_id] = extract.index.to_numpy()
extract = extract.reset_index(drop=True)
extract.to_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(grid_df)
gc.collect()
mean_features = [
'enc_cat_id_mean', 'enc_cat_id_std',
'enc_dept_id_mean', 'enc_dept_id_std',
'enc_item_id_mean', 'enc_item_id_std'
]
df2 = pd.read_feather(f"target_encoding_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")[mean_features]
for store_id in store_id_set_list:
df = pd.read_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
df = pd.concat([df, df2[df2.index.isin(index_store[store_id])].reset_index(drop=True)], axis=1)
df.to_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(df2)
gc.collect()
df3 = pd.read_feather(f"lag_feature_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather").iloc[:, 3:]
for store_id in store_id_set_list:
df = pd.read_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
df = pd.concat([df, df3[df3.index.isin(index_store[store_id])].reset_index(drop=True)], axis=1)
df.to_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
del(df3)
del(store_id_set_list)
gc.collect()
以下函数，相反，只是进一步处理前一个选择，通过删除未使用的特征和重新排序列，并返回用于训练模型的数据：
def load_grid_by_store(end_train_day_x, predict_horizon, store_id):
df = pd.read_feather(f"grid_full_store_{store_id}_{end_train_day_x}_to_{end_train_day_x + predict_horizon}.feather")
remove_features = ['id', 'state_id', 'store_id', 'date', 'wm_yr_wk', 'd', 'sales']
enable_features = [col for col in list(df) if col not in remove_features]
df = df[['id', 'd', 'sales'] + enable_features]
df = reduce_mem_usage(df, verbose=False)
gc.collect()
return df, enable_features
最后，我们现在可以处理训练阶段了。以下代码片段首先定义了训练参数，正如 Monsaraida 所解释的，这是最有效的。由于训练时间的原因，我们只修改了提升类型，选择使用 goss 而不是 gbdt，因为这可以在很大程度上加快训练速度，而不会在性能上损失太多。通过 subsample 参数和特征分数也可以为模型提供良好的加速：在梯度提升的每个学习步骤中，只有一半的示例和一半的特征将被考虑。
此外，在您的机器上使用正确的编译选项编译 LightGBM 也可能提高您的速度，如在这篇有趣的竞赛讨论中所述：
www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/148273
Tweedie 损失，其幂值为 1.1（因此具有与泊松分布相当的基础分布），在建模间歇序列（零销售额占主导地位）时似乎特别有效。所使用的度量标准只是均方根误差（没有必要使用自定义度量标准来表示竞赛度量标准）。我们还使用force_row_wise参数在 Kaggle 笔记本中节省内存。所有其他参数都与 Monsaraida 在其解决方案中提出的参数完全相同（除了 subsampling 参数已被禁用，因为它与 goss 提升类型不兼容）。
Tweedie 损失在哪些 Kaggle 竞赛中已被证明是有用的？你能通过探索 ForumTopics 和 ForumMessages 表来找到关于这种损失及其在 Meta Kaggle 中使用的有用讨论吗？(
www.kaggle.com/datasets/kaggle/meta-kaggle)？
在定义训练参数后，我们只需遍历各个商店，每次上传单个商店的训练数据并训练 LightGBM 模型。每个模型都是经过打包的。我们还从每个模型中提取特征重要性，以便将其合并到一个文件中，然后汇总，从而得到每个特征在该预测范围内的所有商店的平均重要性。
这里是针对特定预测范围训练所有模型的完整函数：
def train(train_df, seed, end_train_day_x, predict_horizon):
lgb_params = {
'boosting_type': 'goss',
'objective': 'tweedie',
'tweedie_variance_power': 1.1,
'metric': 'rmse',
#'subsample': 0.5,
#'subsample_freq': 1,
'learning_rate': 0.03,
'num_leaves': 2 ** 11 - 1,
'min_data_in_leaf': 2 ** 12 - 1,
'feature_fraction': 0.5,
'max_bin': 100,
'boost_from_average': False,
'num_boost_round': 1400,
'verbose': -1,
'num_threads': os.cpu_count(),
'force_row_wise': True,
}
random.seed(seed)
np.random.seed(seed)
os.environ['PYTHONHASHSEED'] = str(seed)
lgb_params['seed'] = seed
store_id_set_list = list(train_df['store_id'].unique())
print(f"training stores: {store_id_set_list}")
feature_importance_all_df = pd.DataFrame()
for store_index, store_id in enumerate(store_id_set_list):
print(f'now training {store_id} store')
grid_df, enable_features = load_grid_by_store(end_train_day_x, predict_horizon, store_id)
train_mask = grid_df['d'] <= end_train_day_x
valid_mask = train_mask & (grid_df['d'] > (end_train_day_x - predict_horizon))
preds_mask = grid_df['d'] > (end_train_day_x - 100)
train_data = lgb.Dataset(grid_df[train_mask][enable_features],
label=grid_df[train_mask]['sales'])
valid_data = lgb.Dataset(grid_df[valid_mask][enable_features],
label=grid_df[valid_mask]['sales'])
# Saving part of the dataset for later predictions
# Removing features that we need to calculate recursively
grid_df = grid_df[preds_mask].reset_index(drop=True)
grid_df.to_feather(f'test_{store_id}_{predict_horizon}.feather')
del(grid_df)
gc.collect()
estimator = lgb.train(lgb_params,
train_data,
valid_sets=[valid_data],
callbacks=[lgb.log_evaluation(period=100, show_stdv=False)],
)
model_name = str(f'lgb_model_{store_id}_{predict_horizon}.bin')
feature_importance_store_df = pd.DataFrame(sorted(zip(enable_features, estimator.feature_importance())),
columns=['feature_name', 'importance'])
feature_importance_store_df = feature_importance_store_df.sort_values('importance', ascending=False)
feature_importance_store_df['store_id'] = store_id
feature_importance_store_df.to_csv(f'feature_importance_{store_id}_{predict_horizon}.csv', index=False)
feature_importance_all_df = pd.concat([feature_importance_all_df, feature_importance_store_df])
pickle.dump(estimator, open(model_name, 'wb'))
del([train_data, valid_data, estimator])
gc.collect()
feature_importance_all_df.to_csv(f'feature_importance_all_{predict_horizon}.csv', index=False)
feature_importance_agg_df = feature_importance_all_df.groupby(
'feature_name')['importance'].agg(['mean', 'std']).reset_index()
feature_importance_agg_df.columns = ['feature_name', 'importance_mean', 'importance_std']
feature_importance_agg_df = feature_importance_agg_df.sort_values('importance_mean', ascending=False)
feature_importance_agg_df.to_csv(f'feature_importance_agg_{predict_horizon}.csv', index=False)
准备好最后一个函数后，我们为我们的管道工作准备好了所有必要的代码。对于封装整个操作的功能，我们需要输入数据集（时间序列数据集、价格数据集、日历信息）以及最后训练日（对于公共排行榜预测为 1,913，对于私有排行榜为 1,941）和预测范围（可能是 7、14、21 或 28 天）。
def train_pipeline(train_df, prices_df, calendar_df,
end_train_day_x_list, prediction_horizon_list):
for end_train_day_x in end_train_day_x_list:
for predict_horizon in prediction_horizon_list:
print(f"end training point day: {end_train_day_x} - prediction horizon: {predict_horizon} days")
# Data preparation
generate_base_grid(train_df, end_train_day_x, predict_horizon)
calc_release_week(prices_df, end_train_day_x, predict_horizon)
generate_grid_price(prices_df, calendar_df, end_train_day_x, predict_horizon)
generate_grid_calendar(calendar_df, end_train_day_x, predict_horizon)
modify_grid_base(end_train_day_x, predict_horizon)
generate_lag_feature(end_train_day_x, predict_horizon)
generate_target_encoding_feature(end_train_day_x, predict_horizon)
assemble_grid_by_store(train_df, end_train_day_x, predict_horizon)
# Modelling
train(train_df, seed, end_train_day_x, predict_horizon)
由于 Kaggle 笔记本有有限的运行时间、有限的内存和磁盘空间，我们建议的策略是复制四个包含此处所示代码的笔记本，并使用不同的预测范围参数进行训练。使用相同的笔记本名称，但部分包含预测参数的值，有助于在另一个笔记本中将模型作为外部数据集收集和处理。
这是第一个笔记本，m5-train-day-1941-horizon-7 (www.kaggle.com/code/lucamassaron/m5-train-day-1941-horizon-7)：
end_train_day_x_list = [1941]
prediction_horizon_list = [7]
seed = 42
train_pipeline(train_df, prices_df, calendar_df, end_train_day_x_list, prediction_horizon_list)
第二个笔记本，m5-train-day-1941-horizon-14 (www.kaggle.com/code/lucamassaron/m5-train-day-1941-horizon-14)：
end_train_day_x_list = [1941]
prediction_horizon_list = [14]
seed = 42
train_pipeline(train_df, prices_df, calendar_df, end_train_day_x_list, prediction_horizon_list)
第三个笔记本，m5-train-day-1941-horizon-21 (www.kaggle.com/code/lucamassaron/m5-train-day-1941-horizon-21)：
end_train_day_x_list = [1941]
prediction_horizon_list = [21]
seed = 42
train_pipeline(train_df, prices_df, calendar_df, end_train_day_x_list, prediction_horizon_list)
最后一个是 m5-train-day-1941-horizon-28 (www.kaggle.com/code/lucamassaron/m5-train-day-1941-horizon-28)：
end_train_day_x_list = [1941]
prediction_horizon_list = [28]
seed = 42
train_pipeline(train_df, prices_df, calendar_df, end_train_day_x_list, prediction_horizon_list)
如果你在一台具有足够磁盘空间和内存资源的本地计算机上工作，你可以一次性运行所有四个预测范围，输入包含它们的列表：[7, 14, 21, 28]。现在，在能够提交我们的预测之前，最后一步是组装预测。
组装公共和私有预测
你可以在这里看到一个关于如何组装公共和私有排行榜预测的示例：
-
公共排行榜示例：
www.kaggle.com/lucamassaron/m5-predict-public-leaderboard -
私人排行榜示例：
www.kaggle.com/code/lucamassaron/m5-predict-private-leaderboard
公共和私人提交之间的变化只是不同的最后训练日：它决定了我们将预测哪些天。
在这个结论性的代码片段中，在加载必要的包，如 LightGBM，对于每个训练结束日和每个预测范围，我们恢复正确的笔记本及其数据。然后，我们遍历所有商店，预测所有商品在从上一个预测范围到现在的范围内的销售额。这样，每个模型都将预测它所训练的单周。
import numpy as np
import pandas as pd
import os
import random
import math
from decimal import Decimal as dec
import datetime
import time
import gc
import lightgbm as lgb
import pickle
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
store_id_set_list = ['CA_1', 'CA_2', 'CA_3', 'CA_4', 'TX_1', 'TX_2', 'TX_3', 'WI_1', 'WI_2', 'WI_3']
end_train_day_x_list = [1913, 1941]
prediction_horizon_list = [7, 14, 21, 28]
pred_v_all_df = list()
for end_train_day_x in end_train_day_x_list:
previous_prediction_horizon = 0
for prediction_horizon in prediction_horizon_list:
notebook_name = f"../input/m5-train-day-{end_train_day_x}-horizon-{prediction_horizon}"
pred_v_df = pd.DataFrame()
for store_index, store_id in enumerate(store_id_set_list):
model_path = str(f'{notebook_name}/lgb_model_{store_id}_{prediction_horizon}.bin')
print(f'loading {model_path}')
estimator = pickle.load(open(model_path, 'rb'))
base_test = pd.read_feather(f"{notebook_name}/test_{store_id}_{prediction_horizon}.feather")
enable_features = [col for col in base_test.columns if col not in ['id', 'd', 'sales']]
for predict_day in range(previous_prediction_horizon + 1, prediction_horizon + 1):
print('[{3} -> {4}] predict {0}/{1} {2} day {5}'.format(
store_index + 1, len(store_id_set_list), store_id,
previous_prediction_horizon + 1, prediction_horizon, predict_day))
mask = base_test['d'] == (end_train_day_x + predict_day)
base_test.loc[mask, 'sales'] = estimator.predict(base_test[mask][enable_features])
temp_v_df = base_test[
(base_test['d'] >= end_train_day_x + previous_prediction_horizon + 1) &
(base_test['d'] < end_train_day_x + prediction_horizon + 1)
][['id', 'd', 'sales']]
if len(pred_v_df)!=0:
pred_v_df = pd.concat([pred_v_df, temp_v_df])
else:
pred_v_df = temp_v_df.copy()
del(temp_v_df)
gc.collect()
previous_prediction_horizon = prediction_horizon
pred_v_all_df.append(pred_v_df)
pred_v_all_df = pd.concat(pred_v_all_df)
当所有预测都收集完毕后，我们使用样本提交文件作为参考将它们合并，包括需要预测的行和列格式（Kaggle 期望在验证或测试期间具有每日销售额的递进列的商品具有不同的行）。
submission = pd.read_csv("../input/m5-forecasting-accuracy/sample_submission.csv")
pred_v_all_df.d = pred_v_all_df.d - end_train_day_x_list
pred_h_all_df = pred_v_all_df.pivot(index='id', columns='d', values='sales')
pred_h_all_df = pred_h_all_df.reset_index()
pred_h_all_df.columns = submission.columns
submission = submission[['id']].merge(pred_h_all_df, on=['id'], how='left').fillna(0)
submission.to_csv("m5_predictions.csv", index=False)
解决方案在私人排行榜上可以达到大约 0.54907，结果是第 12 位，位于金牌区域。恢复 Monsaraida 的 LightGBM 参数（例如，使用 gbdt 而不是 goss 作为提升参数）应该会带来更高的性能（但你需要在本地计算机或 Google Cloud Platform 上运行代码）。
练习
作为练习，尝试比较使用相同迭代次数训练 LightGBM，将提升集设置为 gbdt 而不是 goss。性能和训练时间差异有多大（你可能需要使用本地机器或云机器，因为训练可能超过 12 小时）？
摘要
在本章中，我们面对了一场相当复杂的时间序列竞赛，因此我们尝试的最简单的前几名解决方案实际上相当复杂，它需要编写大量的处理函数。在你阅读完本章后，你应该对如何处理时间序列以及如何使用梯度提升进行预测有一个更好的理解。当数据量足够时，例如这个问题，优先考虑梯度提升解决方案而不是传统方法，应该有助于你为具有层次相关性、间歇序列以及事件、价格或市场条件等协变量可用的问题创建强大的解决方案。在接下来的章节中，你将面对更加复杂的数据竞赛，处理图像和文本。你将惊讶于通过重新创建得分最高的解决方案并理解其内部工作原理，你可以学到多少。
第三章：03 木薯叶病竞赛
加入我们的 Discord 书籍社区
packt.link/EarlyAccessCommunity
在本章中，我们将离开表格数据的领域，专注于图像处理。为了演示在分类竞赛中取得好成绩所需的步骤，我们将使用木薯叶病竞赛的数据：
www.kaggle.com/competitions/cassava-leaf-disease-classification
开始参加 Kaggle 竞赛的第一件事是正确阅读描述：
“作为非洲第二大碳水化合物供应商，木薯是小型农户种植的关键粮食安全作物，因为它能够承受恶劣条件。至少 80%的撒哈拉以南非洲的家庭农场种植这种淀粉根，但病毒性疾病是产量低下的主要原因。借助数据科学，可能能够识别常见疾病，以便进行治疗。”
因此，这个竞赛与一个实际重要的现实生活问题相关——你的里程可能会有所不同，但一般来说，了解这一点是有用的。
“现有的疾病检测方法要求农民寻求政府资助的农业专家的帮助，进行视觉检查和诊断植物。这既费时费力，供应不足，又昂贵。作为一个额外的挑战，有效的解决方案必须在重大限制下表现良好，因为非洲农民可能只能使用移动质量的相机和低带宽。”
这段话——特别是最后一句话——设定了期望：由于数据来自不同的来源，我们可能会遇到一些与图像质量和（可能）分布偏移相关的挑战。
“你的任务是将每个木薯图像分类到四个疾病类别或第五个表示健康叶片的类别。在你的帮助下，农民可能能够快速识别患病植物，从而在它们造成不可修复的损害之前挽救他们的作物。”
这部分相当重要：它指定这是一个分类竞赛，类别数量较少（5 个）。
在完成介绍性的准备工作后，让我们来看看数据。
理解数据和指标
进入这个竞赛的“数据”标签页，我们看到提供的数据摘要：
图 3.1：木薯竞赛数据集描述
我们能从中得到什么？
-
数据格式相当直接，组织者甚至提供了疾病名称和数值代码之间的映射
-
我们有 tfrecord 格式的数据，这对于任何有兴趣使用 TPU 的人来说是个好消息
-
提供的测试集只是一个小的子集，在评估时需要用完整的数据集替换。这表明，在评估时加载先前训练的模型并使用它进行推理是一种更受欢迎的策略。
评估指标被选为分类准确率：developers.google.com/machine-learning/crash-course/classification/accuracy。这个指标接受离散值作为输入，这意味着潜在的集成策略变得稍微复杂一些。损失函数在训练期间实现，以优化学习函数，只要我们想使用基于梯度下降的方法，这个函数就需要是连续的；另一方面，评估指标在训练后用于衡量整体性能，因此可以是离散的。
练习：不构建模型，编写代码进行基本的 EDA
- 比较我们的分类问题中类的基数
通常，这也是检查分布偏移的时刻：如果训练集和测试集中的图像非常不同，这肯定是需要考虑的事情。然而，由于在这种情况下我们没有访问完整的数据集，这一步被省略了——请参阅第一章，其中讨论了对抗性验证，这是一种检测数据集之间概念漂移的流行技术。
构建基线模型
我们从构建基线解决方案开始。运行端到端解决方案的笔记本可在以下位置找到：
虽然希望作为您可能想要尝试的其他比赛的起点是有用的，但遵循本节中描述的流程进行学习更有教育意义，即逐个复制代码单元格，这样您可以更好地理解它（当然，您也可以改进它——它之所以被称为基线解决方案，是有原因的）。
图 3.2：基线解决方案所需的导入
我们首先导入必要的包——虽然个人风格差异是自然的事情，但我们的观点是，将导入集中在一个地方会使代码在比赛进展和向更复杂的解决方案过渡时更容易维护。此外，我们创建了一个配置类：所有定义我们学习过程参数的占位符：
图 3.3：基线解决方案的配置类
组件包括：
-
数据文件夹在您有时在 Kaggle 之外训练模型时非常有用（例如，在 Google Colab 或本地机器上）
-
BATCH_SIZE 是一个参数，如果您想优化您的训练过程（或在受限制的内存环境中处理大图像时使其成为可能），有时需要调整
-
修改 EPOCHS 对于调试很有用：从少量 EPOCHS 开始，以验证你的解决方案是否可以顺畅地从头到尾运行，并在你接近正确解决方案时增加
-
TARGET_SIZE 定义了我们想要将图像重新缩放的尺寸
-
NCLASSES 对应于我们的分类问题中的可能类别数量
编码解决方案的一个好习惯是将重要的部分封装在函数中 - 而创建我们的可训练模型当然符合重要性的标准：
Figure 3.4: 创建我们的模型函数
关于此步骤的一些简要说明：
-
虽然有更多表达性的选项可用，但开始使用一个可以快速迭代的快速模型是实用的；EfficientNet
paperswithcode.com/method/efficientnet架构非常适合这个要求 -
我们添加一个池化层用于正则化目的
-
添加一个分类头 - 一个 Dense 层，其中 CFG.NCLASSES 表示分类器的可能结果数量
-
最后，我们使用与此次比赛要求相对应的损失和指标来编译模型。
练习：检查损失和指标的可能选择 - 一个有用的指南是 keras.io/api/losses/ 其他合理的选项是什么？
下一步是数据：
Figure 3.5: 设置数据生成器
接下来我们设置模型 - 由于我们定义了上述函数，所以非常直接：
Figure 3.6: 实例化模型
在我们开始训练模型之前，我们应该关注回调：
Figure 3.7: 模型回调
值得注意的一些点：
-
ModelCheckpoint 用于确保我们只保留最佳模型的权重，其中最优性是基于要监控的指标（在本例中为验证损失）来决定的
-
EarlyStopping 通过确保在给定数量的 EPOCHS 内验证损失没有改善（下降）时停止训练，帮助我们控制过拟合的风险
-
ReduceLROnPlateau 是一个学习率方案
练习：在上面的设置中，哪些参数值得修改，哪些可以保留为默认值？
使用此设置，我们可以拟合模型：
Figure 3.8: 拟合模型
一旦训练完成，我们可以使用该模型为测试集中的每个图像构建图像类别的预测。回想一下，在这个比赛中，公开（可见）的测试集仅包含一个图像，完整测试集的大小是未知的 - 因此需要稍微复杂的方式来构建提交数据框
Figure 3.9: 生成提交
在本节中，我们展示了如何开始参加一个专注于图像分类的竞赛——你可以使用这种方法快速从基本的 EDA（探索性数据分析）过渡到功能性的提交。然而，像这样的基础方法不太可能产生非常具有竞争力的结果。因此，在下一节中，我们将讨论在顶级得分解决方案中使用的更多专业技术。
顶级解决方案的启示
在本节中，我们汇集了顶级解决方案的各个方面，这些方面可能使我们超越基线解决方案的水平。请注意，在这个竞赛中，排行榜（无论是公开的还是私人的）都非常紧密；这是几个因素的组合：
-
噪声数据——通过正确识别大部分训练数据，很容易达到 0.89 的准确率，然后每个新的正确识别都允许进行微小的提升
-
指标——准确率在集成中可能很棘手
-
数据量有限
预训练
解决数据量有限问题的第一个也是最明显的补救措施是预训练：使用更多数据。木薯竞赛一年前也举行过：
www.kaggle.com/competitions/cassava-disease/overview
通过最小的调整，2019 版的数据可以在当前版本中发挥作用。几位竞争者讨论了这个问题：
-
在论坛上发布了 2019 + 2020 数据集的 TFRecords 格式：
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/199131 -
2019 版竞赛的获胜方案作为一个有用的起点：
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/216985 -
在 2019 年的数据上生成预测并使用伪标签来扩充数据集据报道可以带来一些（轻微的）改进
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/203594
TTA（测试时间增强）
测试时间增强（TTA）背后的想法是对测试图像应用不同的变换：旋转、翻转和平移。这创建了几个不同的测试图像版本，并为每个版本生成一个预测。然后将结果类概率平均，以获得更自信的答案。Andrew Khael 在笔记本中提供了一个关于此技术的优秀演示：www.kaggle.com/code/andrewkh/test-time-augmentation-tta-worth-it
在 Cassava 竞赛中，顶级解决方案广泛使用了 TTA，一个优秀的例子是顶级 3 的私人 LB 结果：www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/221150
Transformer
虽然在竞赛过程中使用了更多的已知架构，如 ResNext 和 EfficientNet，但正是添加了更多新颖的架构，为许多渴望在拥挤的排行榜上取得进步的竞争者提供了额外的优势。Transformer 在 2017 年作为 NLP（自然语言处理）的一次革命性架构出现（如果你错过了启动这一切的论文，这里就是：arxiv.org/abs/1706.03762），并且取得了如此巨大的成功，以至于不可避免地许多人开始思考它们是否也可以应用于其他模态——视觉显然是一个明显的候选者。名为 Vision Transformer 的架构在 Cassava 竞赛中首次亮相。已经公开了一个关于 ViT 的优秀教程：
www.kaggle.com/code/abhinand05/vision-transformer-vit-tutorial-baseline
Vision Transformer 是获胜解决方案的一个重要组成部分：
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/221150
集成学习
集成学习的核心思想在 Kaggle 上非常流行（参见 Kaggle 书籍的第九章，以获取更详细的描述），Cassava 竞赛也不例外。事实证明，结合不同的架构非常有益（通过平均类别概率）：EfficientNet、ResNext 和 ViT 彼此之间足够不同，以至于它们的预测可以相互补充。另一个重要的角度是堆叠，即使用模型分两个阶段：第 3 名解决方案结合了多个模型的预测，然后应用 LightGBM 作为混合器
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/220751
获胜解决方案采用了不同的方法（在最终混合中使用的模型较少），但依赖于相同的核心逻辑：
www.kaggle.com/competitions/cassava-leaf-disease-classification/discussion/221957
摘要
在本章中，我们描述了如何开始使用基线解决方案参加图像分类竞赛，并讨论了多种可能的扩展方法，以进入竞争（获奖）区域。
浙公网安备 33010602011771号