“本文通过 Python 代码实现的方式来来介绍具体的实现。
以 Python 为基础，调用各种比较基础的库，其中穿插一些Kaggle处理的建议，用引号表示。“
1. 数据预处理
1.1 魔术工具及所需要包
魔术工具，python 的 notebook 可以自动的更新 py 文件里的 function。
%load_ext autoreload
%autoreload 2
%matplotlib inline
导入全部的 模块，如果需要安装包的，可以一次性检查一下。
import os
import math
from concurrent.futures import ProcessPoolExecutor
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
1.2 数据导入前的预览
有些竞赛的数据体量巨大，当我们读入 pandas 之后再处理，会大大增加数据处理的时间，因此我们可以选在在 terminal 里面先查阅一下文档的情况。在 notebook 里面使用 ! 可以直接在 cell 里面使用 terminal 的命令行。
!wc -l data/Train.csv
显示文档的长度
---
!shuf -n 1 data/Train.csv
显示文档中的任意一行
---
!head -1 data/Train.csv
显示文档的第 1 行
通过上述预览的方法，我们可以了解文档中各个 columns 的特征之后，来比较好的读取。当我们不告诉 pandas 数据的 types 的时候，他会选择 chunk 的方式来读取，这种方式会占用大量的内存，所以我们尽量指示适合的 types，来加速读取速度，也就是在 read_csv 时，设置 dtypes 参数来帮助 pandas 加速识别。
types = {'id': 'int64',
'item_nbr': 'int32',
'store_nbr': 'int8',
'unit_sales': 'float32',
'onpromotion': 'object'}
%%time df_raw = pd.read_csv(f'{PATH}train.csv', parse_dates=['date'], dtype=types, infer_datetime_format=True)
---
CPU times: user 1min 41s, sys: 5.08s, total: 1min 46s
Wall time: 1min 48s
1.3 nan数据情况及比例预览
list(df_raw.isnull().mean())
---
[0.0, 0.0, 0.0, 0.0, 0.0, 0.050198815830476785, 0.0, 0.6440885010906825, 0.8263907759426613, 0.0, 0.0, 0.0, 0.34201558117793707, 0.8581290121533188, 0.8207067622312246, 0.5254596447491431, 0.0, 0.0, 0.0, 0.0, 0.73982923028981, 0.0008102212527267062, 0.5211542536615769, 0.8027198504206918, 0.6295269554378311, 0.8027198504206918, 0.5432097226550328, 0.8027198504206918, 0.9371293237768775, 0.9371293237768775, 0.9371293237768775, 0.9371293237768775, 0.20082268619507634, 0.9371293237768775, 0.7403876597070739, 0.9371019009037084, 0.9371293237768775, 0.7638691181053288, 0.4666201308818947, 0.8916597070738548, 0.8918990339669679, 0.8918990339669679, 0.7528127142411967, 0.7510202555313181, 0.7526506699906513, 0.7524761607977563, 0.7526506699906513, 0.7528127142411967, 0.8038716110937987, 0.8009772514802119, 0.8009747584917419, 0.826959177313805, 0.8270638828295419]
解释：isnull 函数计算 True/False，之后用 mean 函数来计算 True/（True+False）的数量，也就是 nan 的占比。
df_raw.describe(include = 'all').T
---
describe函数可以查阅所有特征的情况，是一种比较好的概览方式。
1.4 对 datetime 处理特征
def add_datetime(df, col):
ts = ['time','year','month','day','hour','minute','second','week','weekofyear','dayofweek', 'weekday','dayofyear','quarter','is_month_start','is_month_end','is_quarter_start', 'is_quarter_end','is_year_start','is_year_end','daysinmonth']
for item in ts: df[item] = getattr(df[col].dt, item)
return df.drop(col, axis=1)
%prun add_datetime(df_raw, 'saledate')
由于时间序列文档里面，重要的特征包含了日期、周中第几天、月中第几天、第几周、是否是周末等数据。因此我们就用这个编写的 add_datetime 函数来加入这些特征。
1.5 处理 category
def create_cat(df):
cols = df.select_dtypes(include='object').columns
for col in cols: df[col] = df[col].astype('category')
%time create_cat(df_raw)
为 category 建立编码，但是这里的编码创建之后，是将 column 的类型改称 'category'，还不是机器学习模型可以处理的数字编码。
def create_num(df):
df_cat = df.copy()
cols = df.select_dtypes(include='category')
for col in cols:
df_cat[col] = df_cat[col].cat.codes
return df_cat
%time df_cat = create_num(df_raw)
因此，通过 create_num 函数来进行编码的数字化，但是仍然需要保留原始的 df_raw，再最后解释模型的时候，这些编码的原始对应关系需要被反映设回到原来的样子。
1.6 存储格式
简单的格式转化完后，我们可以先存一下数据，往往我们会为了 pandas 的读取方便，将文件存储成 feather 格式。
对 category 的映射文档，我们存成 csv 格式，另一个可以用于机器学习的文档，我们存成 feather 格式。
os.makedirs('tmp', exist_ok=True)
cols = df_raw.select_dtypes(include='category')
cols.to_csv('tmp/df_raw')
df_cat.to_feather('tmp/df_cat')
2. 模型实现
2.1 数据划分成训练集和验证集
之前我们介绍过，为了模拟时间序列模型的预测未来特征，我们按照时间顺序把训练集和验证集切开，这样验证集的评判就是对未来的预测。
df_nan = df_cat.dropna()
%time x, y = sample_split(df_nan, 'saledate', 'SalePrice', 2000)
%time x_train, x_val, y_train, y_val = train_val_split(x, y, 1600)
去掉模型不可以处理的 nan 函数，然后切分训练集和验证集。
*其实这里直接去掉nan并不是很好，可能会丢失掉一些特征，我们之后再来优化。
2.2 评价函数实现
def rmse(x, y):
return math.sqrt(((x-y)**2).mean())
def print_score(m):
res = { 'rmse_train': rmse(m.predict(x_train), y_train),
'rmse_val':rmse(m.predict(x_val), y_val),
'score_train':m.score(x_train, y_train),
'score_val': m.score(x_val, y_val)}
if hasattr(m, 'oob_score_'): res['oob_score'] = m.oob_score_
return res
这里的 rmse 是 kaggle 模型中的要求，所以我们自己实现一下。
2.3 模型运行及调试
model = RandomForestRegressor(n_estimators=20, oob_score=True)
model.fit(x_train, y_train)
print_score(model)
---
{'rmse_train': 0.13393435504176915, 'rmse_val': 0.34025337669453626, 'score_train': 0.9644767465836186, 'score_val': 0.7752985602391357, 'oob_score': 0.749322023225568}
继续优化一下：
model = RandomForestRegressor(n_estimators=40, min_samples_leaf=3, max_features=0.5, oob_score=True)
model.fit(x_train, y_train)
print_score(model)
---
{'rmse_train': 0.18659299774817495, 'rmse_val': 0.33602501198638357, 'score_train': 0.9310523550318314, 'score_val': 0.7808486360579396, 'oob_score': 0.7597670944592627}
hmm，结果看上去一般般，下次我们来研究一下如何优化。
3. 模型可靠性评估
3.1 模型的可靠性
由于 RF 是很多的树，也就是 estimators 很多，类似模型集合的构造。因此，当我们获取每个 estimators 的估算，并评价其 std 标准方差，我们就能够获得可靠性。
%time preds = np.stack([t.predict(x_val) for t in model.estimators_])
print(np.mean(preds), np.std(preds))
---
CPU times: user 63 ms, sys: 0 ns, total: 63 ms
Wall time: 62.4 ms
10.173505185913063 0.6872819393762499
那这样，我们就能得到更多对模型精度的认知了。
*3.2 提高速度的并行计算
在 RF 中，我们所有的树都是独立的，基于此假设我们很容易得到，每棵树在预测期都是可以独立运算的。因此我们在这里可以采用并行计算的方法来实现，通过 python 的多线程管理包。
from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor
def parallel_trees(m, func, n_jobs=8):
return list(ProcessPoolExecutor(n_jobs).map(func, m.estimators_))
def parallel_preds(t):
return t.predict(x_val)
---
%time preds = np.stack(parallel_trees(model, parallel_preds))
print(np.mean(preds), np.std(preds))
查看一下 preds 的 shape，我们会发现这是一个（40，400）的 matrix，表示了每个 tree 对都 400 个 validation 数据提供了预测，因此通过这些数据我们可以来进一步针对不同的模型特征进行可靠性评估。
3.3 模型不同特征的可靠性评估，基于不同的组别
如果我们需要评估不同特征在不同组别上的可靠性，那我们就可以进一步的通过 pandas 自带的 groupby 功能来进行合理的模型精度上的分析。
Eg，我们的模型里面有一个特征是 UsageBand，我们就来使用这个特征来看看，该模型在不同组别上的精度情况。
sum_val = x_val.copy()
sum_val['pre_mean'] = np.mean(preds, axis=0)
sum_val['pre_std'] = np.std(preds, axis=0)
sum_val.UsageBand.value_counts().plot.barh() para = ['UsageBand', 'pre_mean', 'pre_std'] sum_val_ = sum_val[para].groupby('UsageBand', as_index=False).mean() sum_val_.plot('UsageBand', 'pre_mean', 'barh', xerr='pre_std', alpha=0.6, xlim=(8,11))
(sum_val_.pre_mean/sum_val_.pre_std).sort_values(ascending=False)
同样，也可以通过 sort_values 的方式来输出根据 groupby 之后，可信度比较低的类别，作为对于模型使用及解释的认知。
4. 特征重要性 Feature Importance
RF 中特别重要的一个参数是 Feature Importance，这个参数可以告诉我们在所有的 model feature 里面，哪些对于最终的结果比较重要，类似与 PCA 的分析结果。
def rf_feat_importance(m, df):
return pd.DataFrame({'cols':df.columns, 'imp':m.feature_importances_} ).sort_values('imp', ascending=False)
---
fi = rf_feat_importance(model, x_train); fi[:10]
这个工具最大的作用，就是帮你理解你的模型特征本身是不是存在未来函数，这个解释是不是合理，让你的机器学习变得透明。当然许多 Kaggle 竞赛就是在这些数据里面出现了 data leakage 的情况的。但是一般来说，选手都会公开这些特征，让大家再次站在统一起跑线上，所以要一直跟踪你的 Kaggle forum。
此外，去掉冗余的数据也不会加速我们模型的训练，比如我们只选取 imf 高于 0.005 的特征来进行 RF 的训练。
解释一下这个贡献：
RF 对于贡献的定义就是，针对这个 feature，我们 shuffle 一下，看看最终的 score 降低了多少，就是这个 feature 的贡献。
