【实战篇】随机森林预测气温（二）

回顾

我们先来回顾 【实战篇】随机森林预测气温（一） 中完成的内容。小鱼读取了 2016 年全年的气温数据集，包含 8 和特征和 1 个标签，样本个数也比较少：

>> df.shape
(348, 9)

之后，为了观察特征的情况，我们将特征中的 year month day 组合成日期类型，通过绘制折线图观察了 temp_1 temp_2 friend 以及真实气温 actual 随日期的变化情况。

接下来，小鱼进行了简单的数据预处理，将 week 属性类型的特征转换为度热编码的形式。

最后将数据集按照 1:3 划分为测试集和训练集，使用随机森林建模，并借助绝对百分比误差 MAPE 评估该回归预测任务。

>> print(f'MAPE:{mape.mean():.2%}')
MAPE:6.08%

并且使用可视化的方式，将测试集的预测结果和真实气温的分布直观地呈现出来。

本节，我们主要讨论增大数据量、引入新特征对结果的影响。

读取新的数据集

读取新文件 temps_extended.csv ：

import pandas as pd
import os

df = pd.read_csv("data" + os.sep + "temps_extended.csv")
df.head()

新的数据中，数据规模发生了变化，由 348 增加到了 2191 ，范围从 2016 年变为 2011-2016 年。

>> df.shape
(2191, 12)
>> df.year.value_counts().sort_index()
2011    365
2012    365
2013    365
2014    365
2015    365
2016    365
2017      1
Name: year, dtype: int64

并且加入了新的天气指标：

• ws_1 前一天的风速
• prcp_1 前一天的降水
• snwd_1 前一天的积雪深度

既然有了新的特征，先来看看他们长什么样吧~

组合日期：

# 处理时间数据
from datetime import datetime

# datetime 格式
dates = [datetime(year,month,day) 
         for year,month,day in zip(df.year, df.month, df.day)]
dates[:5]

绘制特征和日期的折线图：

import matplotlib.pyplot as plt

%matplotlib inline

def feature_plot(ax, feature, xlable=''):
    ax.plot(dates, df[feature], linewidth=0.6)
    ax.set_xlabel(xlable)
    ax.set_title(feature)

plt.style.use("seaborn-whitegrid")
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(10,16), dpi=80)
fig.autofmt_xdate(rotation = 45)

# 标签值、昨天、前天、噪声、平均最高气温、风速、降水、积雪
for ax, feature in zip(axes.flatten(), ('actual','temp_1','temp_2', 'friend','average','ws_1','prcp_1','snwd_1')):
    feature_plot(ax, feature)

从时间上来看，前一天的风速 ws_1 和噪音很相似，是非常不稳定的。前一天有积雪 snwd_1 的天数屈指可数。

合成新的特征

在数据分析和特征提取的过程中，我们的出发点都是尽可能多的选择有价值的特征，因为我们能得到的信息越多，之后建模可以利用的信息也是越多的。

比如在这份数据中，我们有完整的日期数据，但是显示天气的变换肯定是跟季节因素有关的。可原始数据集中并没有体现出季节的指标，我们可以自己创建一个季节变量当做新的特征，无论是对之后建模还是分析都会起到帮助的。

def check_season(month):
    if month in [1,2,12]:
        return 'winter'
    elif month in [3,4,5]:
        return 'spring'
    elif month in [6,7,8]:
        return 'summer'
    elif month in [9,10,11]:
        return 'fall'
    
df['season'] = df.month.apply(check_season)
df

添加季节后的数据集：

各季节包含的记录数：

>> df.season.value_counts()
spring    552
summer    552
fall      546
winter    541
Name: season, dtype: int64

有了季节特征之后，假如想观察一下不同季节的时候上述各项指标的变换情况该怎么做呢？

这里给大家推荐一个非常实用的绘图函数 pairplot ，需要我们先安装 seaborn 这个工具包，它相当于是在 Matplotlib 的基础上进行封装，说白了就是用起来更简单规范了。

import seaborn as sns

sns.set(style="ticks", color_codes=True)

sns.pairplot(
    df[['temp_1', 'prcp_1', 'average', 'actual', 'season']], 
    hue="season", 
    diag_kind="kde", 
    palette=sns.color_palette("hls", 4), 
    plot_kws=dict(alpha=0.7), 
    diag_kws=dict(shade=True)
)

PairPlot 绘制结果：

image.png

PairPlot 绘图中，x 轴和 y 轴都是我们的 4 个特征，不同颜色的点表示不同的季节。在主对角线上 x 轴和 y 轴都是同一个特征，因此绘制的是该特征的核密度估计图，表示当前特征在不同季节时的数值分布情况。

注：核密度估计图可以理解为直方图中，将每个柱子分得非常细，然后用一根平滑的曲线把无穷个小柱子的顶端用线连起来，是一个微分再积分的过程。

PairPlot 绘图的其他位置则用散点图来表示两个特征之间的关系。例如在左下角 temp_1 和 actual 就呈现出了很强的正相关性，相关系可以用 y=x 这样的直线来拟合；prcp_1 和 autual 也呈现出一定的相关性，但相关性较弱。

扩大数据集、增加特征的影响

值增大数据集，特征保持和原始小的数据集一致，结果会有提升吗？

定义相关函数，实现数据集切分、随机森林模型训练以及预测：

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor

def get_train_test(df):
    data = pd.get_dummies(df)
    y = data.actual
    X = data.drop('actual',axis=1)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=0)
    return X_train, X_test, y_train, y_test

def fit_and_predict(x_train, y_train, x_test, y_test, print_mape=True):
    
    rfr = RandomForestRegressor(n_estimators=100, random_state=0)
    rfr.fit(x_train, y_train)
    predictions = rfr.predict(x_test)
    # 计算误差 MAPE (Mean Absolute Percentage Error)
    errors = abs(predictions - y_test)
    mape = errors / y_test
    if print_mape:
        print(f'MAPE:{mape.mean():.2%}')
    return rfr, mape.mean()

读取原始数据集：

small_df = pd.read_csv('data/temps.csv')
small_df = pd.get_dummies(small_df, prefix={'week':'weekday'})

切分新的数据集以及原始小数据集：

>> X_train, X_test, y_train, y_test = get_train_test(df)
>> X_train.shape, X_test.shape
((1643, 21), (548, 21))
>> X_train_s, X_test_s, y_train_s, y_test_s = get_train_test(small_df)
>> X_train_s.shape, X_test_s.shape
((261, 14), (87, 14))

为了测试效果能够公平，使用同一个测试集 X_test：

>> small_df_features = X_train_s.columns.array
>> tree_s, mape_s = fit_and_predict(X_train_s, y_train_s, X_test[small_df_features], y_test)
MAPE:7.93%
>> tree_s_feature, mape_s_feature = fit_and_predict(X_train[small_df_features], y_train, X_test[small_df_features], y_test)
MAPE:6.88%

可以看到，当我们把数据量增大之后，绝对百分比误差从 7.93% 下降到了 6.88%。在机器学习任务中，我们希望数据量能够越大越好，这样可利用的信息就更多了，机器学习任务数据量一般在几万到几百万。

下面再对比一下特征数量对结果的影响，之前这两次比较还没有加入新的特征，这回我们把降水，风速，积雪和我们自己构造的季节特征加入训练集中，看看效果又会怎样：

>> tree, mape = fit_and_predict(X_train, y_train, X_test, y_test)
MAPE:6.63%

模型整体效果有了略微提升，误差从 6.88% 下降到 6.63%，下降了 0.25 个百分点。

特征重要性分析

这回特征也多了，我们可以好好研究下特征重要性这个指标了。

获取特征及特征重要性，并进行就地降序排列：

feature_importance = pd.Series(
    tree.feature_importances_,
    index=tree.feature_names_in_
)
feature_importance.sort_values(ascending=False, inplace=True)
feature_importance

特征及特征重要性：

temp_1           0.848102
average          0.050617
ws_1             0.020563
friend           0.016345
temp_2           0.014183
day              0.013133
year             0.009602
prcp_1           0.007477
month            0.004651
weekday_Fri      0.002502
weekday_Sun      0.001914
weekday_Tues     0.001842
weekday_Wed      0.001776
weekday_Mon      0.001753
weekday_Sat      0.001699
weekday_Thurs    0.001168
season_summer    0.000878
season_spring    0.000756
season_fall      0.000702
season_winter    0.000238
snwd_1           0.000100
dtype: float64

对结果进行可视化展示：

feature_importance.plot(kind='bar', color='r', edgecolor='k', linewidth=1.2)
plt.title('Features Importances')

特征重要性如何加以利用呢？

答案就特征的累加重要性：先把特征按照其重要性进行排序，再算起累计值，通常我们都以 95% 为阈值，看看有多少个特征累加在一起之后，其特征重要性的累加值超过该阈值，就取它们当做筛选后的特征。

注：计算累加使用 cumsum() 函数：比如 cumsum([1,2,3,4]) 得到的结果就是其累加值 (1,3,6,10)。

>> importance_cumsum = feature_importance.cumsum()
>> importance_cumsum
temp_1           0.848102
average          0.898720
ws_1             0.919282
friend           0.935627
temp_2           0.949810
day              0.962943
year             0.972545
prcp_1           0.980023
month            0.984674
weekday_Fri      0.987176
weekday_Sun      0.989089
weekday_Tues     0.990932
weekday_Wed      0.992707
weekday_Mon      0.994461
weekday_Sat      0.996159
weekday_Thurs    0.997327
season_summer    0.998204
season_spring    0.998960
season_fall      0.999663
season_winter    0.999900
snwd_1           1.000000
dtype: float64

可视化展示特征重要性累加值：

importance_cumsum.plot(color='r', linewidth=1.2)

# 在 0.95 的位置画一条红色虚线
plt.hlines(
    y=0.95,
    xmin=0,
    xmax=len(importance_cumsum),
    linestyles='dashed',
    color='g',
)

plt.xticks(
    np.linspace(0, len(importance_cumsum)-1, len(importance_cumsum)), 
    importance_cumsum.index, 
    rotation='vertical'
)

plt.title('Cumulative Importances')

这里当第 6 个特征出现的时候，其总体的累加值超过了 95%，达到96.29%。那么接下来我们的对比实验又来了，如果只用这 6 个特征来计算 MAPE。

>> feature_6 = np.array(importance_cumsum[:6].index)
>> feature_6
array(['temp_1', 'average', 'ws_1', 'friend', 'temp_2', 'day'],
      dtype=object)
>> tree_feature, mape_feature = fit_and_predict(X_train[feature_6], y_train, X_test[feature_6], y_test)
MAPE:6.78%

使用累加重要性超过 95% 的特征建模，损失从 6.63% 增加到了 6.78% ，误差上升了 0.15 个百分点。

注：随机森林的算法本身就会考虑特征的问题，会优先选择有价值的，我们认为的去掉一些，相当于可供候选的就少了，出现这样的现象在随机森林中并不奇怪！

时间效率分析

虽然模型没有提升，我们还可以再看看在时间效率的层面上有没有进步。定义函数 calculate_time_mape 计算 10 次训练、预测预测过程的平均耗时。

import time

def calculate_time_mape(features, x_train, y_train, x_test, y_test):
    times = []
    for _ in range(10):
        start_time = time.time()
        rfr, mape = fit_and_predict(x_train[features], y_train, x_test[features], y_test, print_mape=False)
        end_time = time.time()
        times.append(end_time - start_time)
    return round(mape,4), round(np.mean(times),2)

分别计算原始数据集、增加数据量之后的数据集、增加数据量之后并减少特征的数据集耗时情况：

model_comparison = pd.DataFrame(
    columns=['mape','run_time'], 
    index=['original','exp_all','exp_reduced']
)

model_comparison.loc['original'] = calculate_time_mape(small_df_features, X_train_s, y_train_s, X_test,y_test)
model_comparison.loc['exp_all'] = calculate_time_mape(X_train.columns.array, X_train, y_train, X_test, y_test)
model_comparison.loc['exp_reduced'] = calculate_time_mape(feature_6, X_train, y_train, X_test, y_test)
model_comparison

耗时、MAPE 情况：

其中：

• original 代表老数据，也就是量少特征少的那份；
• exp_all 代表完整新数据；
• exp_reduced 代表按照 95% 阈值选择的部分重要特征数据集。

将对比结果进行可视化展示：

fontdict = {'fontsize': 18}
fontdict_yaxis = {'fontsize': 14}

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,5), sharex=True)
plt.xticks(np.linspace(0,2,3), model_comparison.index)
x_values = model_comparison.index

# 误差对比
ax1.bar(x_values, model_comparison.mape*100, color=['b', 'r', 'g'], edgecolor='k', linewidth=1.5)
ax1.set_ylabel('MAPE (%)', fontdict=fontdict_yaxis)
ax1.set_title('Model Error Comparison', fontdict=fontdict)

# 时间效率对比
ax2.bar(x_values, model_comparison.run_time, color=['b', 'r', 'g'], edgecolor='k', linewidth=1.5)
ax2.set_ylabel('Run Time (sec)', fontdict=fontdict_yaxis) 
ax2.set_title('Model Run-Time Comparison', fontdict=fontdict)

绘图结果：

从柱形图来看，选取重要特征建模虽然误差有小幅上升，但在时间上却下降明显，牺牲了小的准确率，换来了效率上的提升。