一、前言
回归随机森林作为一种机器学习和数据分析领域常用且有效的算法，对其原理和代码实现过程的掌握是非常有必要的。为此，本文将着重介绍从零开始实现回归随机森林的过程，对于随机森林和决策树的相关理论原理将不做太深入的描述。本文的目的只是为了演示回归随机森林主要功能的具体实现过程，在实现过程中不会考虑代码性能，会更加注重代码可读性。
实现语言：Python
依赖：pandas, numpy
二、原理介绍
随机森林属于Bagging类算法，而Bagging 又属于集成学习一种方法（集成学习方法大致分为Boosting和Bagging方法，两个方法的不同请参考[10]），集成学习的大致思路是训练多个弱模型打包起来组成一个强模型，强模型的性能要比单个弱模型好很多（三个臭皮匠顶一个诸葛亮。注意：这里的弱和强是相对的），其中的弱模型可以是决策树、SVM等模型，在随机森林中，弱模型选用决策树。
在训练阶段，随机森林使用bootstrap采样从输入训练数据集中采集多个不同的子训练数据集来依次训练多个不同决策树；在预测阶段，随机森林将内部多个决策树的预测结果取平均得到最终的结果。本文主要介绍回归随机森林从零实现的过程，实现的RFR(回归随机森林)具有以下功能：
- 模型训练
- 模型数据预测
- 计算feature importance
2.1 模型训练
2.1.1 基础理论
本文实现的RFR是将多个二叉决策树（即CART，这也是sklearn,spark内部实现的模型）打包组合而成的，训练RFR便是训练多个二叉决策树。在训练二叉决策树模型的时候需要考虑怎样选择切分变量(特征）、切分点以及怎样衡量一个切分变量、切分点的好坏。针对于切分变量和切分点的选择，本实现采用穷举法，即遍历每个特征和每个特征的所有取值，最后从中找出最好的切分变量和切分点；针对于切分变量和切分点的好坏，一般以切分后节点的不纯度来衡量，即各个子节点不纯度的加权和
其中，
- 备注1：
Gini不纯度只适用于分类任务，
- 备注2：
信息熵只适用于分类任务， 其他同备注1。
- 备注3
MSE只适用于回归任务，
- 备注4
MAE只适用于回归任务，
- 备注5
在sklearn内部，DecisionTreeClassifier, RandomForestClassifier等基于决策树的分类模型默认使用'gini'作为impurity function，也可通过criterion参数指定为'entropy' ；而DecisionTreeRegressor, RandomForestRegressor等基于决策树的回归模型默认使用'mse'作为impurity function，也可通过criterion参数指定为'mae'。
决策树中某一节点的训练过程在数学等价于下面优化问题：
即寻找
在本文实现的回归随机森林中，
2.1.2 训练过程
RFR训练过程示意图如图2所示。
说明：
- bootstrap[1]是对输入训练样本集合 N次放回重复抽样得到样本集合
进行
N一般等于(
sklearn的实现代码)，这暗示了的大小[2]，具体可以查看
中的某些样本会在
中重复出现多次。在对
进行bootstrap采样时，
63.3%的概率被抽中。为什么是63.3%？考虑N次放回重复抽样，每次抽样每个样本被选中的概率为中每个样本会以
，进行N次抽样，被选中概率为：
具体可以参考[3-4]。
- subsample会根据输入参数sample_sz的大小从 sample_sz个样本组成subsample样本集合
中采样
sklearn中，subsample大小就是输入训练样本集合(在
的大小N，不可以通过输入参数指定subsample的大小，详见这里)。
- 为什么要进行bootstrap? 集成学习（包括随机森林）的目的是为了使用多个不同的子模型来增加最终模型预测结果的鲁棒性和稳定性(即减小方差)，如果多个子模型都采用同样的数据集训练，那么训练得出的子模型都是相同的，集成学习将变得没有意义，为了能从原输入训练数据集得到多个不同的数据集，需要使用bootstrap来从原数据集采样出不同的子数据集来训练子模型。
图2中n个回归决策树的训练过程如图3所示。
2.2 模型预测
RFR的预测结果是由内部所有二叉决策树的预测结果取平均值得到的。二叉决策树的预测过程主要分为以下步骤：
- 针对某一输入样本，从二叉决策树的根节点起，判断当前节点是否为叶子节点，如果是则返回叶子节点的预测值(即当前叶子中样本目标变量的平均值），如果不是则进入下一步；
- 根据当前节点的切分变量的和切分值，将样本中对应变量的值与节点的切分值对比。如果样本变量值小于等于当前节点切分值，则访问当前节点的左子节点；如果样本变量值大于当前节点切分值，则访问当前节点的右子节点；
- 循环步骤2，直到访问到叶子节点，并返回叶子节点的预测值。
2.3 计算feature importance
特征的重要性表示特征对预测结果影响程度，某一特征重要性越大，表明该特征对预测结果的影响越大，重要性越小，表明该特征对预测结果越小。RFR中某一特征的重要性，是该特征在内部所有决策树重要性的平均值，而在决策树中，计算某一个特征的重要性可以采用以下三种方法：
方法1
步骤如下：
- 使用训练数据训练模型；
- 计算训练数据在模型上依据某一metric的score，记为
（在回归中，可以选用r2）；
- 遍历训练数据集中的每一个feature，每次在原训练数据集的基础上将对应的feature进行shuffle操作，然后使用模型得到shuffle后数据集的score，记为
，最后通过
计算出第
个feature的重要性。
方法2
如果一个feature很重要，那么从数据集中去除该feature后，模型的prediction error会增大；相反，如果一个feature不重要，那么从数据集中去除后，模型的prediction error不会变化太大。该方法的步骤如下：
- 使用原数据集训练模型，并评估模型在训练数据集上的prediction error，记为
;
- 遍历训练数据集中的每一个feature，每次从原训练数据集的基础上去除该feature，然后使用得到数据集训练模型，最终计算出prediction error，记为
，最后通过
计算出第
个feature的重要性。
该方法和方法1有点类似，但该方法在计算每个feature的重要性的时候都需要重新训练模型，会增加额外的计算量，在实际应用一般不采用该方法。
方法3
该方法是sklearn内部树模型计算feature重要性时采用的方法。即某一节点
其中，
为了使所有feature的重要性加起来等于1，需要每一feature的重要性进行normalization，即公式(2-7)
方法3在sklearn中的实现，请查看_tree.pyx。
在本文实现的RFR中，同时实现了方法1和方法3。
三、回归随机森林实现
3.1 代码
代码有点长，不想看的可以直接跳过。
from multiprocessing import Pool, cpu_countclass RFRegressor(object):"""Random Forest Regressor"""def __init__(self, n_trees, sample_sz, min_leaf_sz=5, n_jobs=None, max_depth=None):self._n_trees = n_treesself._sample_sz = sample_szself._min_leaf_sz = min_leaf_szself._n_jobs = n_jobsself._max_depth = max_depthself._trees = [self._create_tree() for i in range(self._n_trees)]def _get_sample_data(self, bootstrap=True):"""Generate training data for each underlying decision treeParameters----------bootstrap: boolean value, True/FalseThe default value is True, it would bootstrap sample from input training data. If False, the exclusive sampling willbe performed.Returns-------idxs: array-like objectReturn the indices of sampled data from input training data"""if bootstrap:idxs = np.random.choice(len(self._X), self._sample_sz)else:idxs = np.random.permutation(len(self._X))[:self._sample_sz]return idxsdef _create_tree(self):"""Build decision treeeReturns-------dtree : DTreeRegressor object"""return DTreeRegressor(self._min_leaf_sz, self._max_depth)def _single_tree_fit(self, tree):"""Fit the single underlying decision treeParameters----------tree : DTreeRegressor objectReturns-------tree : DTreeRegressor object"""sample_idxs = self._get_sample_data()return tree.fit(self._X.iloc[sample_idxs], self._y[sample_idxs])def fit(self, x, y):"""Train a forest regressor of trees from the training set(x, y)Parameters----------x : DataFrame,The training input samples.y : Series or array-like objectThe target values."""self._X = xself._y = yif self._n_jobs:self._trees = self._parallel(self._trees, self._single_tree_fit, self._n_jobs)else:for tree in self._trees:self._single_tree_fit(tree)def predict(self, x):"""Predict target values using trained modelParameters---------x : DataFrame or array-like objectinput samplesReturns-------ypreds : array-like objectpredicted target values"""all_tree_preds = np.stack([tree.predict(x) for tree in self._trees]) return np.mean(all_tree_preds, axis=0)def _parallel(self, trees, fn, n_jobs=1):"""Parallel jobs executionParameters----------trees : list-like objecta list-like object contains all underlying treesfn : function-like objectmap functionn_jobs : integerThe number of jobs.Returns-------result : list-like objecta list-like result object for each call of map function `fn`"""try:workers = cpu_count()except NotImplementedError:workers = 1if n_jobs:workders = n_jobspool = Pool(processes=workers)result = pool.map(fn, trees)return list(result)@propertydef feature_importances_(self):"""Calculate feature importanceReturns-------self._feature_importances : array-like objectthe importance score of each feature """if not hasattr(self, '_feature_importances'):norm_imp = np.zeros(len(self._X.columns))for tree in self._trees:t_imp = tree.calc_feature_importance()norm_imp = norm_imp + t_imp / np.sum(t_imp)self._feature_importances = norm_imp / self._n_treesreturn self._feature_importances@propertydef feature_importances_extra(self):"""Another method to calculate feature importance"""norm_imp = np.zeros(len(self._X.columns))for tree in self._trees:t_imp = tree.calc_feature_importance_extra()norm_imp = norm_imp + t_imp / np.sum(t_imp)norm_imp = norm_imp / self._n_treesimp = pd.DataFrame({'col':self._X.columns, 'imp':norm_imp}).sort_values('imp', ascending=False)return impclass DTreeRegressor(object):def __init__(self, min_leaf_sz, max_depth=None):self._min_leaf_sz = min_leaf_szself._split_point = 0self._split_col_idx = 0self._score = float('inf')self._sample_sz = 0self._left_child_tree = Noneself._right_child_tree = Noneself._feature_importances = []self._node_importance= 0if max_depth is not None:max_depth = max_depth - 1self._max_depth = max_depthdef fit(self, x, y):self._X = xself._y = yself._col_names = self._X.columnsself._feature_importances = np.zeros(len(self._col_names))self._sample_sz = len(self._X)self._val = np.mean(self._y)if self._max_depth is not None and self._max_depth < 2:return selfself._find_best_split()return selfdef _calc_mse_inpurity(self, y_squared_sum, y_sum, n_y):"""Calculate Mean Squared Error impurityThis is just the recursive version for calculating varianceParameters----------y_squared_sum: float or int , the sum of y squared y_sum: float or int , the sum of y valuen_y: int, the number of samplesReturns-------"""dev = (y_squared_sum / n_y) - (y_sum / n_y) ** 2return devdef _find_best_split(self):for col_idx in range(len(self._col_names)):self._find_col_best_split_point(col_idx)self._feature_importances[self._split_col_idx] = self._node_importanceif self.is_leaf:return left_child_sample_idxs = np.nonzero(self.split_col <= self.split_point)[0]right_child_sample_idxs = np.nonzero(self.split_col > self.split_point)[0]self._left_child_tree = (DTreeRegressor(self._min_leaf_sz, self._max_depth).fit(self._X.iloc[left_child_sample_idxs], self._y[left_child_sample_idxs]))self._right_child_tree = (DTreeRegressor(self._min_leaf_sz, self._max_depth).fit(self._X.iloc[right_child_sample_idxs], self._y[right_child_sample_idxs]))def _find_col_best_split_point(self, col_idx):x_col = self._X.values[:, col_idx]sorted_idxs = np.argsort(x_col)sorted_x_col = x_col[sorted_idxs]sorted_y = self._y[sorted_idxs]lchild_n_samples = 0lchild_y_sum = 0.0lchild_y_square_sum = 0.0rchild_n_samples = self._sample_szrchild_y_sum = sorted_y.sum()rchild_y_square_sum = (sorted_y ** 2).sum()node_y_sum = rchild_y_sumnode_y_square_sum = rchild_y_square_sum for i in range(0, self._sample_sz - self._min_leaf_sz):xi, yi = sorted_x_col[i], sorted_y[i]rchild_n_samples -= 1rchild_y_sum -= yirchild_y_square_sum -= (yi ** 2)lchild_n_samples += 1lchild_y_sum += yilchild_y_square_sum += (yi ** 2)if i < self._min_leaf_sz or xi == sorted_x_col[i+1]:continuelchild_impurity = self._calc_mse_inpurity(lchild_y_square_sum,lchild_y_sum, lchild_n_samples)rchild_impurity = self._calc_mse_inpurity(rchild_y_square_sum,rchild_y_sum, rchild_n_samples)split_score = (lchild_n_samples * lchild_impurity + rchild_n_samples * rchild_impurity) / self._sample_szif split_score < self._score:self._score = split_scoreself._split_point = xiself._split_col_idx = col_idxself._node_importance = (self._sample_sz * (self._calc_mse_inpurity(node_y_square_sum, node_y_sum, self._sample_sz) - split_score))def predict(self, x):if type(x) == pd.DataFrame:x = x.valuesreturn np.array([self._predict_row(row) for row in x])def _predict_row(self, row):if self.is_leaf:return self._valt = (self._left_child_tree if row[self._split_col_idx] <= self.split_point else self._right_child_tree)return t._predict_row(row)def __repr__(self):pr = f'sample: {self._sample_sz}, value: {self._val}rn'if not self.is_leaf:pr += f'split column: {self.split_name}, split point: {self.split_point}, score: {self._score} 'return prdef calc_feature_importance(self):if self.is_leaf:return self._feature_importancesreturn (self._feature_importances + self._left_child_tree.calc_feature_importance()+ self._right_child_tree.calc_feature_importance())def calc_feature_importance_extra(self):imp = []o_preds = self.predict(self._X.values)o_r2 = metrics.r2_score(self._y, o_preds)for col in self._col_names:tmp_x = self._X.copy()shuffle_col = tmp_x[col].valuesnp.random.shuffle(shuffle_col)tmp_x.loc[:,col] = shuffle_coltmp_preds = self.predict(tmp_x.values)tmp_r2 = metrics.r2_score(self._y, tmp_preds)imp.append((o_r2 - tmp_r2))imp = imp / np.sum(imp)return imp@propertydef split_name(self):return self._col_names[self._split_col_idx]@propertydef split_col(self):return self._X.iloc[:, self._split_col_idx]@propertydef is_leaf(self):return self._score == float('inf')@propertydef split_point(self):return self._split_point3.2 测试
训练与预测：
使用sklearn内部的RandomForestRegressor的结果：
需要注意的是，上面两次训练的样本不同，所以和实际情况有点出入，但大体趋势是对。
feature importance:
RFRegressor采用方法3计算的feature重要性：
RandomForestRegressor计算的feature重要性：
RFRegressor采用方法1计算的feature重要性：
参考与荐读：
[1] Bootstrap Sample: Definition, Example
[2] Number of Samples per-Tree in a Random Forest
[3] Why on average does each bootstrap sample contain roughly two thirds of observations?
[4] Random Forests out-of-bag sample size
[5] Random Forests Leo Breiman and Adele Cutler
[6] Why is Random Forest with a single tree much better than a Decision Tree classifier?
[7] Winning the bias-variance tradeoff
[8] subsample-bootstrapping
[9] How can SciKit-Learn Random Forest sub sample size may be equal to original training data size?
[10] the difference between bagging and boosting
[11] The Mathematics of Decision Trees, Random Forest and Feature Importance in Scikit-learn and Spark
[12] Sklearn Decision Trees documentation
