第9章--随机森林项目实战——气温预测（1/2）
第8章已经讲解过随机森林的基本原理，本章将从实战的角度出发，借助Python工具包完成气温预测任务，其中涉及多个模块，主要包含随机森林建模、特征选择、效率对比、参数调优等。这个例子实在太长了，分为三篇介绍。这是第一篇。
随机森林建模：气温预测的任务目标就是使用一份天气相关数据来预测某一天的最高温度，属于回归任务，首先观察一下数据集：
输出结果中表头的含义如下。
- year,moth,day,week：分别表示的具体的时间。
- temp_2：前天的最高温度值。
- temp_1：昨天的最高温度值。
- average：在历史中，每年这一天的平均最高温度值。
- actual：就是标签值，当天的真实最高温度。
- friend：这一列可能是凑热闹的，你的朋友猜测的可能值，不管它就好。
该项目实战主要完成以下3项任务。
1.使用随机森林算法完成基本建模任务：包括数据预处理、特征展示、完成建模并进行可视化展示分析。
2.分析数据样本量与特征个数对结果的影响：在保证算法一致的前提下，增加数据样本个数，观察结果变化。重新考虑特征工程，引入新特征后，观察结果走势。
3.对随机森林算法进行调参，找到最合适的参数：掌握机器学习中两种经典调参方法，对当前模型选择最合适的参数。
9.1.1特征可视化与预处理
拿到数据之后，一般都会看看数据的规模，做到心中有数：
print('数据维度:', features.shape) #数据维度: (348, 9)
输出结果显示该数据一共有348条记录，每个样本有9个特征。如果想进一步观察各个指标的统计特性，可以用.describe()展示：
输出结果展示了各个列的数量，如果有数据缺失，数量就会有所减少。由于各列的统计数量值都是348，所以表明数据集中并不存在缺失值，并且均值、标准差、最大值、最小值等指标都在这里显示。
对于时间数据，也可以进行格式转换，原因在于有些工具包在绘图或者计算的过程中，用标准时间格式更方便：
1 # 处理时间数据 2 import datetime 3 4 # 分别得到年，月，日 5 years = features['year'] 6 months = features['month'] 7 days = features['day'] 8 9 # datetime格式 10 dates = [str(int(year)) + '-' + str(int(month)) + '-' + str(int(day)) for year, month, day in zip(years, months, days)] 11 dates = [datetime.datetime.strptime(date, '%Y-%m-%d') for date in dates]
为了更直观地观察数据，最简单有效的办法就是画图展示，首先导入Matplotlib工具包，再选择一个合适的风格（其实风格差异并不是很大）：
1 # 准备画图 2 import matplotlib.pyplot as plt 3 4 %matplotlib inline 5 6 # 指定默认风格 7 plt.style.use('fivethirtyeight')
开始布局，需要展示4项指标，分别为最高气温的标签值、前天、昨天、朋友预测的气温最高值。既然是4个图，不妨采用2×2的规模，这样会更清晰，对每个图指定好其图题和坐标轴即可：
1 # 设置布局 2 fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize = (10,10)) 3 fig.autofmt_xdate(rotation = 45) 4 5 # 标签值 6 ax1.plot(dates, features['actual']) 7 ax1.set_xlabel(''); ax1.set_ylabel('Temperature'); ax1.set_title('Max Temp') 8 9 # 昨天 10 ax2.plot(dates, features['temp_1']) 11 ax2.set_xlabel(''); ax2.set_ylabel('Temperature'); ax2.set_title('Previous Max Temp') 12 13 # 前天 14 ax3.plot(dates, features['temp_2']) 15 ax3.set_xlabel('Date'); ax3.set_ylabel('Temperature'); ax3.set_title('Two Days Prior Max Temp') 16 17 # 我的逗逼朋友 18 ax4.plot(dates, features['friend']) 19 ax4.set_xlabel('Date'); ax4.set_ylabel('Temperature'); ax4.set_title('Friend Estimate') 20 21 plt.tight_layout(pad=2)
上述代码可以生成图9-1的输出。
图9-1 各项特征指标
由图可见，各项指标看起来还算正常（由于是国外的天气数据，在统计标准上有些区别）。接下来，考虑数据预处理的问题，原始数据中的week列并不是一些数值特征，而是表示星期几的字符串，计算机并不认识这些数据，需要转换一下。
图9-2是常用的转换方式，称作one-hot encoding或者独热编码，目的就是将属性值转换成数值。对应的特征中有几个可选属性值，就构造几列新的特征，并将其中符合的位置标记为1，其他位置标记为0。
图9-2 特征编码
既可以用Sklearn工具包中现成的方法完成转换，也可以用Pandas中的函数，综合对比后觉得用Pandas中的.get_dummies()函数最容易：
1 # 独热编码 2 features = pd.get_dummies(features) 3 features.head(5)
完成数据集中属性值的预处理工作后，默认会把所有属性值都转换成独热编码的格式，并且自动添加后缀，这样看起来更清晰。
其实也可以按照自己的方式设置编码特征的名字，在使用时，如果遇到一个不太熟悉的函数，想看一下其中的细节，一个更直接的方法，就是在Notebook中直接调用help工具来看一下它的API文档，下面返回的就是get_dummies的细节介绍，也可以查阅在线文档：
help(pd.get_dummies) Help on function get_dummies in module pandas.core.reshape.reshape: get_dummies(data, prefix=None, prefix_sep='_', dummy_na=False, columns=None, sparse=False, drop_first=False, dtype=None) -> 'DataFrame' Convert categorical variable into dummy/indicator variables. Parameters ---------- data : array-like, Series, or DataFrame Data of which to get dummy indicators. prefix : str, list of str, or dict of str, default None String to append DataFrame column names. Pass a list with length equal to the number of columns when calling get_dummies on a DataFrame. Alternatively, `prefix` can be a dictionary mapping column names to prefixes. prefix_sep : str, default '_' If appending prefix, separator/delimiter to use. Or pass a list or dictionary as with `prefix`. dummy_na : bool, default False Add a column to indicate NaNs, if False NaNs are ignored. columns : list-like, default None Column names in the DataFrame to be encoded. If `columns` is None then all the columns with `object` or `category` dtype will be converted. sparse : bool, default False Whether the dummy-encoded columns should be backed by a :class:`SparseArray` (True) or a regular NumPy array (False). drop_first : bool, default False Whether to get k-1 dummies out of k categorical levels by removing the first level. dtype : dtype, default np.uint8 Data type for new columns. Only a single dtype is allowed. .. versionadded:: 0.23.0 Returns ------- DataFrame Dummy-coded data. See Also -------- Series.str.get_dummies : Convert Series to dummy codes. Examples -------- >>> s = pd.Series(list('abca')) >>> pd.get_dummies(s) a b c 0 1 0 0 1 0 1 0 2 0 0 1 3 1 0 0 >>> s1 = ['a', 'b', np.nan] >>> pd.get_dummies(s1) a b 0 1 0 1 0 1 2 0 0 >>> pd.get_dummies(s1, dummy_na=True) a b NaN 0 1 0 0 1 0 1 0 2 0 0 1 >>> df = pd.DataFrame({'A': ['a', 'b', 'a'], 'B': ['b', 'a', 'c'], ... 'C': [1, 2, 3]}) >>> pd.get_dummies(df, prefix=['col1', 'col2']) C col1_a col1_b col2_a col2_b col2_c 0 1 1 0 0 1 0 1 2 0 1 1 0 0 2 3 1 0 0 0 1 >>> pd.get_dummies(pd.Series(list('abcaa'))) a b c 0 1 0 0 1 0 1 0 2 0 0 1 3 1 0 0 4 1 0 0 >>> pd.get_dummies(pd.Series(list('abcaa')), drop_first=True) b c 0 0 0 1 1 0 2 0 1 3 0 0 4 0 0 >>> pd.get_dummies(pd.Series(list('abc')), dtype=float) a b c 0 1.0 0.0 0.0 1 0.0 1.0 0.0 2 0.0 0.0 1.0
特征预处理完成之后，还要把数据重新组合一下，特征是特征，标签是标签，分别在原始数据集中提取一下：
print('Shape of features after one-hot encoding:', features.shape) #Shape of features after one-hot encoding: (348, 15)
1 # 数据与标签 2 import numpy as np 3 4 # 标签 5 labels = np.array(features['actual']) 6 7 # 在特征中去掉标签 8 features= features.drop('actual', axis = 1) 9 10 # 名字单独保存一下，以备后患 11 feature_list = list(features.columns) 12 13 # 转换成合适的格式 14 features = np.array(features)
在训练模型之前，需要先对数据集进行切分：
1 # 数据集切分 2 from sklearn.model_selection import train_test_split 3 4 train_features, test_features, train_labels, test_labels = train_test_split(features, labels, test_size = 0.25, 5 random_state = 42) 6 print('训练集特征:', train_features.shape) 7 print('训练集标签:', train_labels.shape) 8 print('测试集特征:', test_features.shape) 9 print('测试集标签:', test_labels.shape)
训练集特征: (261, 14) 训练集标签: (261,) 测试集特征: (87, 14) 测试集标签: (87,)
9.1.2随机森林回归模型
万事俱备，开始建立随机森林模型，首先导入工具包，先建立1000棵树模型试试，其他参数暂用默认值，然后深入调参任务：
1 # 导入算法 2 from sklearn.ensemble import RandomForestRegressor 3 4 # 建模 5 rf = RandomForestRegressor(n_estimators= 1000, random_state=42) 6 7 # 训练 8 rf.fit(train_features, train_labels)
RandomForestRegressor(bootstrap=True, ccp_alpha=0.0, criterion='mse', max_depth=None, max_features='auto', max_leaf_nodes=None, max_samples=None, min_impurity_decrease=0.0, min_impurity_split=None, min_samples_leaf=1, min_samples_split=2, min_weight_fraction_leaf=0.0, n_estimators=1000, n_jobs=None, oob_score=False, random_state=42, verbose=0, warm_start=False)
由于数据样本量非常小，所以很快可以得到结果，这里选择先用MAPE指标进行评估，也就是平均绝对百分误差。
1 # 预测结果 2 predictions = rf.predict(test_features) 3 4 # 计算误差 5 errors = abs(predictions - test_labels) 6 7 # mean absolute percentage error (MAPE) 8 mape = 100 * (errors / test_labels) 9 10 print ('MAPE:',np.mean(mape))
MAPE: 6.011244187972058
其实对于回归任务，评估方法还是比较多的，下面列出几种，都很容易实现，也可以选择其他指标进行评估。
9.1.3树模型可视化方法
得到随机森林模型后，现在介绍怎么利用工具包对树模型进行可视化展示，首先需要安装Graphviz工具，其配置过程如下。
第① 步：下载安装。
登录网站https://graphviz.gitlab.io/_pages/Download/Download_windows.html
下载graphviz-2.38.msi，完成后双击这个msi文件，然后一直单击next按钮，即可安装Graphviz软件（注意：一定要记住安装路径，因为后面配置环境变量会用到路径信息，系统默认的安装路径是C:Program Files (x86)Graphviz2.38）。
第②步：配置环境变量。
将Graphviz安装目录下的bin文件夹添加到Path环境变量中。
本例中：
D: oolsGraphVizin
第③步：验证安装。
进入Windows命令行界面，输入“dot–version”命令，然后按住Enter键，如果显示Graphviz的相关版本信息，则说明安装配置成功，
dot -version dot - graphviz version 2.38.0 (20140413.2041) libdir = "D: oolsGraphVizin" Activated plugin library: gvplugin_dot_layout.dll Using layout: dot:dot_layout Activated plugin library: gvplugin_core.dll Using render: dot:core Using device: dot:dot:core The plugin configuration file: D: oolsGraphVizinconfig6 was successfully loaded. render : cairo dot fig gd gdiplus map pic pov ps svg tk vml vrml xdot layout : circo dot fdp neato nop nop1 nop2 osage patchwork sfdp twopi textlayout : textlayout device : bmp canon cmap cmapx cmapx_np dot emf emfplus eps fig gd gd2 gif gv imap imap_np ismap jpe jpeg jpg
metafile pdf pic plain plain-ext png pov ps ps2 svg svgz tif tiff tk vml vmlz vrml wbmp xdot xdot1.2 xdot1.4 loadimage : (lib) bmp eps gd gd2 gif jpe jpeg jpg png ps svg
最后还需安装graphviz、pydot和pydotplus插件，在命令行中输入相关命令即可，代码如下：
1 pip3 install graphviz 2 pip3 install pydot2 3 pip3 install pydotplus 4 pip3 install pydot
上述工具包安装完成之后，就可以绘制决策树模型：
1 # 导入所需工具包 2 from sklearn.tree import export_graphviz 3 import pydot #pip install pydot 4 5 # 拿到其中的一棵树 6 tree = rf.estimators_[5] 7 8 # 导出成dot文件 9 export_graphviz(tree, out_file = 'tree.dot', feature_names = feature_list, rounded = True, precision = 1) 10 11 # 绘图 12 (graph, ) = pydot.graph_from_dot_file('tree.dot') 13 14 # 展示 15 graph.write_png('tree.png');
执行完上述代码，会在指定的目录下（如果只指定其名字，会在代码所在路径下）生成一个tree.png文件，这就是绘制好的一棵树的模型，如图9-8所示。树模型看起来有点太大，观察起来不太方便，可以使用参数限制决策树的规模，还记得剪枝策略吗？预剪枝方案在这里可以派上用场。
1 print('The depth of this tree is:', tree.tree_.max_depth) 2 #The depth of this tree is: 15
图9-9对生成的树模型中各项指标的含义进行了标识，看起来还是比较好理解，其中非叶子节点中包括4项指标：所选特征与切分点、评估结果、此节点样本数量、节点预测结果（回归中就是平均）。
图9-9 树模型可视化中各项指标含义
9.1.4特征重要性
讲解随机森林算法的时候，曾提到使用集成算法很容易得到其特征重要性，在sklearn工具包中也有现成的函数，调用起来非常容易：
1 # 得到特征重要性 2 importances = list(rf.feature_importances_) 3 4 # 转换格式 5 feature_importances = [(feature, round(importance, 2)) for feature, importance in zip(feature_list, importances)] 6 7 # 排序 8 feature_importances = sorted(feature_importances, key = lambda x: x[1], reverse = True) 9 10 # 对应进行打印 11 [print('Variable: {:20} Importance: {}'.format(*pair)) for pair in feature_importances]
Variable: temp_1 Importance: 0.7 Variable: average Importance: 0.19 Variable: day Importance: 0.03 Variable: temp_2 Importance: 0.02 Variable: friend Importance: 0.02 Variable: month Importance: 0.01 Variable: year Importance: 0.0 Variable: week_Fri Importance: 0.0 Variable: week_Mon Importance: 0.0 Variable: week_Sat Importance: 0.0 Variable: week_Sun Importance: 0.0 Variable: week_Thurs Importance: 0.0 Variable: week_Tues Importance: 0.0 Variable: week_Wed Importance: 0.0
上述输出结果分别打印了当前特征及其所对应的特征重要性，绘制成图表分析起来更容易：
1 # 转换成list格式 2 x_values = list(range(len(importances))) 3 4 # 绘图 5 plt.bar(x_values, importances, orientation = 'vertical') 6 7 # x轴名字 8 plt.xticks(x_values, feature_list, rotation='vertical') 9 10 # 图名 11 plt.ylabel('Importance'); plt.xlabel('Variable'); plt.title('Variable Importances');
上述代码可以生成图9-10的输出，可以明显发现，temp_1和average这两个特征的重要性占据总体的绝大部分，其他特征的重要性看起来微乎其微。那么，只用最厉害的特征来建模，其效果会不会更好呢？其实并不能保证效果一定更好，但是速度肯定更快，先来看一下结果：
图9-10 随机森林特征重要性
1 # 选择最重要的那两个特征来试一试 2 rf_most_important = RandomForestRegressor(n_estimators= 1000, random_state=42) 3 4 # 拿到这俩特征 5 important_indices = [feature_list.index('temp_1'), feature_list.index('average')] 6 train_important = train_features[:, important_indices] 7 test_important = test_features[:, important_indices] 8 9 # 重新训练模型 10 rf_most_important.fit(train_important, train_labels) 11 12 # 预测结果 13 predictions = rf_most_important.predict(test_important) 14 15 errors = abs(predictions - test_labels) 16 17 # 评估结果 18 19 mape = np.mean(100 * (errors / test_labels)) 20 21 print('mape:', mape)
mape: 6.229055723613811
从损失值上观察，并没有下降，反而上升了，说明其他特征还是有价值的，不能只凭特征重要性就否定部分特征数据，一切还要通过实验进行判断。
但是，当考虑时间效率的时候，就要好好斟酌一下是否应该剔除掉那些用处不大的特征以加快构建模型的速度。到目前为止，已经得到基本的随机森林模型，并可以进行预测，下面来看看模型的预测值与真实值之间的差异：
1 # 日期数据 2 months = features[:, feature_list.index('month')] 3 days = features[:, feature_list.index('day')] 4 years = features[:, feature_list.index('year')] 5 6 # 转换日期格式 7 dates = [str(int(year)) + '-' + str(int(month)) + '-' + str(int(day)) for year, month, day in zip(years, months, days)] 8 dates = [datetime.datetime.strptime(date, '%Y-%m-%d') for date in dates] 9 10 # 创建一个表格来存日期和其对应的标签数值 11 true_data = pd.DataFrame(data = {'date': dates, 'actual': labels}) 12 13 # 同理，再创建一个来存日期和其对应的模型预测值 14 months = test_features[:, feature_list.index('month')] 15 days = test_features[:, feature_list.index('day')] 16 years = test_features[:, feature_list.index('year')] 17 18 test_dates = [str(int(year)) + '-' + str(int(month)) + '-' + str(int(day)) for year, month, day in zip(years, months, days)] 19 20 test_dates = [datetime.datetime.strptime(date, '%Y-%m-%d') for date in test_dates] 21 22 predictions_data = pd.DataFrame(data = {'date': test_dates, 'prediction': predictions}) 23 24 # 真实值 25 plt.plot(true_data['date'], true_data['actual'], 'b-', label = 'actual') 26 27 # 预测值 28 plt.plot(predictions_data['date'], predictions_data['prediction'], 'ro', label = 'prediction') 29 plt.xticks(rotation = '60'); 30 plt.legend() 31 32 # 图名 33 plt.xlabel('Date'); plt.ylabel('Maximum Temperature (F)'); plt.title('Actual and Predicted Values');
通过上述输出结果的走势可以看出，模型已经基本能够掌握天气变化情况，接下来还需要深入数据，考虑以下几个问题。
1.如果可利用的数据量增大，会对结果产生什么影响呢？
2.加入新的特征会改进模型效果吗？此时的时间效率又会怎样？
未完待续。
该书资源下载，请至异步社区：https://www.epubit.com