Python代码实践|随机森林是“黑盒子”？不存在的，撸完代码你就懂了

点击上方“AI公园”，关注公众号

---

作者：Sonali Dasgupta

编译：ronghuaiyang

动机: 随机森林集成的广泛应用在实际的问题中，既有分类也有回归。它们的热门可以归功于使用随机森林方法使用者只需要很少的数据清洗，不需要进行特征的缩放，就可以得到最优的结果。为了更好的理解背后的原理，我决定从头实现一个随机森林，然后在数据的不同的特征上进行一些可视化，看看这些特征在决定最终结果中扮演了什么样的角色。人们常常认为随机森林是个黑盒子，但是简单的过一遍这个算法，就会发现，这除了是一个利用“多数投票”方法的强大的技术外，实际上可解释性还是挺强的。

介绍

随机森林集成是一个分而治之的方法，用来提高单个弱决策树模型的能力。背后的主要原则是通过一组“弱学习器”，组合起来形成一个“强学习器”。每个分类器都是一个独立的“弱学习器”，所有的分类器一起组成一个“强学习器”。

随机森林集成

随机森林的核心概念就是对一些决策树的结果取平均。随机森林是对那些不熟悉机器学习的人解释预测背后的直觉的一个常用的工具。但是解释随机森林如何得到这个预测，使用哪些特征或者独立变量，还是个麻烦事。随机森林常常被错误的解释为一个“黑盒子”，很难理解。

Random Forest Ensemble

用Python写一个随机森林

下面是一个从头写随机森林代码和优化的练习，可以增进对随机森林的深入理解。我们还会看看，如何理解给定的特征的权值是如何给到的，为什么这么设置，为什么不用其他的特征来预测结果。

随机森林和决策树类

在Python中使用面向对象编程，定义一个类RandomForest，包含了DecisionTree的对象。随机森林基本上就是做一个收集和平局的计算，每一个核心方程使用决策树来构造。

class RandomForest:	
    def __init__(self, sample_size, n_trees = 10 , min_samples_leaf = 5, max_depth = 4, max_features = None):	
        self.n_trees, self.sample_size , self.min_samples_leaf, self.max_depth, self.max_features = n_trees, sample_size, min_samples_leaf, max_depth, max_features	
        self.trees = [self.create_tree() for _ in range(self.n_trees)]	
    def create_tree(self):	
        return DecisionTree(min_samples_leaf = self.min_samples_leaf, max_depth = self.max_depth, max_features = self.max_features)	
    def fit(self, X, y):   	
        #calls the var_split method of underlying Decision Trees	
        for tree in self.trees:	
            random_idxs = np.random.permutation(X.shape[0])[:self.sample_size]	
            tree.var_split(X.iloc[random_idxs, :], y[random_idxs]) 	
    def predict(self, x):	
      #  average of the predictions from each tree	
        return np.mean([t.predict(x) for t in self.trees], axis = 0)	
    def plot_pdp(self, X, y, feature_name, n_clusters = 0): pass	
    def plot_pdp(self, X, y, feature_names, n_clusters = 0): pass	
    def find_feature_importances(self, X, y): pass

The Random Forest class

class DecisionTree:	
    def __init__(self, int min_samples_leaf = 3, int max_depth = 4, int level = 0, max_features = None, float parent_value = float('inf')):	
        #initialize Decision Tree	
    @property 	
    def is_leaf(self):	
        #a node is a leaf if its score is equal to infinity	
    def find_better_split(self, x, y, int ind): 	
        """logic to find the best value to split the tree 	
        on for a particular feature"""	
    def var_split(self, x, y):	
      """logic to find the best feature to split the tree on, and	
      its optimum value for the split"""	
    def predict_row(self, row):      	
        """logic to predict the result for 	
        a single data sample (row) in the	
        test / validation set"""	
    def predict_row_for_ti(self, row, feat_contribs):	
        #functionality for Tree Interpreter	
    def predict(self, X):	
        y_pred = []	
        for row in range(X.shape[0]):	
            y_pred.append(self.predict_row(X.iloc[row, :]))	
        return y_pred	
    def get_prediction_and_bias(self):	
         #returns the prediction of this tree and bias ( the value at the root )	
    def get_child_trees(self):	
        return self.leftTree, self.rightTree	
    def __repr__(self):	
        return "score: " +str(self.score) + " avg: "+str(self.value) +  " split val: " + str(self.split_val) + " split feature : "+ str(self.split_feature)tTree.predict_row_for_ti(row, feat_contribs)

The Decision Tree class

• 构造器（init），输入参数包括每棵树的样本数量，每棵树分到的最大特征数量，目的是为每棵树注入随机性，移除预测的偏差。输入参数还包括一些随机森林算法用到的其他的超参数。

• 拟合函数实现了关键的功能，使用训练数据构建随机森林中的每棵树。

• 预测方法实现了对测试数据进行每一棵树的预测和平均的功能，得到测试数据的最终的结果。

训练模型的逻辑流程

Random Forest类中的fit()方法用来训练模型，调用了 Decision Tree 类中的varsplit()函数，又在循环中调用了 findbetter_split()。这样构造了随机森林中的每一个决策树成员，对每棵树使用训练数据进行拟合。

def var_split(self, x, y):	
        if x.shape[0] > self.min_samples_leaf and self.level < self.max_depth - 1:	
            #keep traversing till conditions for a leaf node are met	
            if self.max_features is not None:	
                if self.max_features in ['auto', 'sqrt']:	
                    self.max_features = int(math.sqrt(x.shape[1]))	
                else:	
                    if self.max_features == 'log2':	
                        self.max_features = int(np.log(float(x.shape[1]))/np.log(2))	
                    else:	
                        if isinstance(self.max_features, float):	
                            self.max_features = int(self.max_features * x.shape[1])	
                        else:	
                            self.max_features = x.shape[1]	
            else:	
                self.max_features = x.shape[1]	
            self.max_features = int(self.max_features)	
            feature_inds = np.random.permutation(x.shape[1])[:self.max_features] 	
            #iterate over all features	
            feature_inds = [index for index in feature_inds if x.columns[index] != None]	
            for ind in feature_inds:	
                self.find_better_split(x, y, ind) 	
           #finds best feature to split on 	
            if self.parent_value == float('inf'):	
                self.parent_value  = self.value	
            x_lhs, x_rhs = x.iloc[:self.split_val,:], x.iloc[self.split_val:,:]	
            self.leftTree = DecisionTree(min_samples_leaf = self.min_samples_leaf, max_depth = self.max_depth, level = self.level + 1, parent_value = self.parent_value)	
            self.leftTree.var_split(x_lhs, y[:self.split_val])	
            self.rightTree = DecisionTree(min_samples_leaf = self.min_samples_leaf, max_depth = self.max_depth, level = self.level + 1, parent_value = self.parent_value)	
            self.rightTree.var_split(x_rhs, y[self.split_val:])	
        else :	
            self.score = float('inf')	
            #samples are averaged to form final prediction at leaf node	
            y = [val for val in y if val != None]	
            self.value = np.mean(y)

var_split() function implementation in Decision Tree

 def find_better_split(self, x, y, int ind):	
        x1 = x.values[:, ind]	
        inds = np.argsort(x1, axis = 0)	
        sorted_y, sorted_x = y[inds], x1[inds]	
        rhs_count, rhs_sum, rhs_sum2  = x.shape[0], sorted_y.sum(), (sorted_y**2).sum() 	
        lhs_count, lhs_sum, lhs_sum2 = 0, 0., 0.	
      #  find optimal value to split of this feature the tree on	
        for i in range(0, x.shape[0]-self.min_samples_leaf+1):	
            lhs_count, lhs_sum, lhs_sum2 = lhs_count+1, lhs_sum + sorted_y[i], lhs_sum2 + sorted_y[i]**2	
            rhs_count, rhs_sum, rhs_sum2 = rhs_count-1, rhs_sum - sorted_y[i], rhs_sum2 - sorted_y[i]**2	
            if i < self.min_samples_leaf - 1 or sorted_x[i] == sorted_x[i+1]:	
                continue	
            updated_score = ((lhs_count * stddev(lhs_count, lhs_sum2, lhs_sum)) + (rhs_count * stddev(rhs_count, rhs_sum2, rhs_sum)))/(x.shape[0])	
            #updated score is the sum of standard deviations in left and right subtrees, which need to be minimized	
            if updated_score < self.score :	
                self.score = updated_score	
                self.split_feature = x.columns[ind]	
                self.split_val = i	
                self.value = (np.mean(y[:i])*i + np.mean(y[i:])*(x.shape[0] - i))/(x.shape[0])	
        self.score = self.score

findbettersplit() function implementation

通过寻找分割之后最大化同质性的特征，决策树递归的分成左右两个子树，这通过两个循环完成：

• 第一个循环遍历所有的特征，决定那个是最好的进行分割。函数 varsplit()对每一个特征调用 findbetter_split() ，看看这个特征是否能够得到最纯的叶子节点。

• findbettersplit() 中的内部的循环中，为每一个特征计算最优化的分割值。

• 计算一个得分，也就是在不同的值上分割每个特征的信息增益。这个得分表示了分割之后的子树的多样性，需要最小化。

• 选出来进行分割的特征和特征的值的组合会让子树有最低的多样性，或者说最大的同质性。

• 通过寻找最佳的特征和最佳的分割值，树分成了左右两个子树，左树和右树的类的成员增加了。

• 在左树和右树成员上递归的调用 var_split()函数，找到最好的特征和最佳的分割值，继续下去，直到叶子节点的其中一个条件满足，比如说最大的深度或者是这个树已经长完了，或者数据样本没有了。

• 这个叶子节点的值就是所包含的所有的训练样本的均值，或者是这叶子节点所对应的值。

对测试/验证数据集做预测

 def predict(self, x):	
        #averages the predictions of each of its Decision Trees	
        return np.mean([t.predict(x) for t in self.trees], axis = 0)

predict() function for Random Forest

def predict(self, X):	
        #append the predictions for each data sample or row by calling predict_row()	
        y_pred = []	
        for row in range(X.shape[0]):	
            y_pred.append(self.predict_row(X.iloc[row, :]))	
        return y_pred	
def predict_row(self, row):      	
        if self.is_leaf: 	
            #return the value at the leaf node as a prediction for this row	
            return self.value	
        if row[self.split_feature] < self.split_val:	
            return self.leftTree.predict_row(row)	
        else:	
            return self.rightTree.predict_row(row)

predict() implementation for a Decision Tree

预测结果的代码非常好理解：

• 随机森林的集成预测是通过对每一个独立的决策树求均值得到的。

注意：随机森林回归器的平均值预测是个很好的技术，但是对于随机森林分类器，考虑到树的预测的频率，通常选择最大频率的预测值。

• 在决策树内部，预测每一个样本，通过遍历树，对比特征的值，对每个level进行分割，然后走到左树或者右树，直到到达叶子节点。

• 对样本的预测存在叶子节点中，在训练的时候已经提前被赋值了。

解释随机森林模型

我们现在有了随机森林集成的基本的结构。我们再加一些非常有用的增强功能，特征交互和树的解释，这些对于解释随机森林是非常重要的概念。

特征重要性背后的直觉

import numpy as np	
def find_feature_importances(self, X, y):	
    feat_imp = {}	
    y_pred_old = self.predict(X)	
    for ind in range(X.shape[1]):	
        X_new = X.copy()	
        np.random.shuffle(X_new.values[:, ind])	
        y_pred = self.predict(X_new)	
        feat_imp[X.columns[ind]] = np.fabs(np.sum(y_pred_old) - np.sum(y_pred))	
        del X_new	
    """Keeping all other feature values constant, shuffle this feature and 	
    calculate the difference between the earlier results and the results using 	
    random values for this feature.	
    A high difference denotes high feature importance."""	
    features, importances =zip(*sorted(feat_imp.items(), reverse=True))	
    return features, importances	
RandomForest.find_feature_importances = find_feature_importances	
idxs = np.random.permutation(X_train.shape[0])[:50]	
features, imp = find_feature_importances(rf_mine, X_train.iloc[idxs, :], y_train[idxs])	
plt.barh(features[:10], imp[:10])

Calculate feature importances for this model

为了找到一个特定的特征对于模型的预测有多少决定性，下面是一个方法：

• 在数据集上的预测的评估度量（准确率或者误差）记录下来。

• 接下来，开始决定特征的值的重要性，使用 shuffle() 进行随机的打乱。

• 使用shuffle之后的这个特定的特征，重新运行这个模型，对数据集进行预测，重新评估。

• 计算评估的结果和没有打乱之前的结果的差。

• 如果差别很大，结果变得很差，这个特定的特征就对模型的预测具有很大的决定性。

如果一个特征的随机打乱的值给模型带来了巨大的性能下降，我们说这个用来分割的特征在树中的位置相当的高，而且重要。在我的实现中，每个特征都会使用这种技术。

可以用条状图来可视化特征以及对于的重要性，用来解释我们的随机森林模型的不同的特征的权值。

Feature Importance Chart

随机森林的局部依赖绘图

局部依赖绘图对于理解特征和目标之间的关系特别的有用。是单调的上升，下降，还是不规则的？

单变量绘图 vs 局部依赖绘图

如果我们看一个典型的单变量的图，去找特征对于目标变量的作用，我们很少会发现连续的上升/下降的曲线或者直线，单变量的图往往是不规则的。

举个例子，比如你有一个房价数据集，预测房价，给你的是一组特征，包括建房年月，地点等其它特征。如果你画一个单变量的图，建房日期和，你期待的是一个连续的上升的曲线，因为房价随着建房日期应该是上升的。但是画出来的图却是在某些年之间有下降。这个可能会错误的引导你认为这些年建的房子对于房价来说是不好的年。

Univariate Plot for the log of House Prices varying with Year Built

然而，可能是正好那些年买的多数的房子在一个相对便宜的地区，导致了曲线的下降。

局部依赖绘图会解决这个问题，展示房价和建房时间的真正关系，你会发现确实是持续上升的。这要归功于局部依赖绘图，我们让考虑的特征中只有一个变化，让其他的特征保持不变。这样我们去掉了其他的依赖变量的噪声，得到了一个对于每个特征和目标变量的真实的自然的解释。

Partial Dependence Plot would show the correct trend of prices against Year Built

局部依赖绘图的实现

我对于局部依赖绘图的实现如下：

import numpy as np	
import matplotlib.pyplot as plt	
def plot_pdp(self, X, y, feature_name, n_clusters = 2):	
    feature_id = -1	
    y_pred_vals = []	
    y_pred_mean = []	
    for index in range(X.shape[1]):	
        if X.columns[index] == feature_name:	
            feature_id = index	
    X_new = X.copy()	
    for feature_val in np.unique(X.values[:, feature_id]):	
        X_new.values[:, feature_id] = feature_val	
        y_pred_val = self.predict(X_new)	
        y_pred_vals.append(y_pred_val)   	
    del X_new    	
    """keeping the values of all other features constant, find the value 	
    of target variable for each uniqye value of this feature"""	
    y_pred_randoms = np.random.normal(y_pred_vals)[:n_clusters]       	
    y_pred_mean = np.mean(y_pred_vals, axis = 0)	
    print("y_pred_mean: "+str(y_pred_mean))	
    plt.plot(X.values[:, feature_id], y_pred_mean)	
    plt.show()	
RandomForest.plot_pdp = plot_pdp	
idxs = np.random.permutation(X_train.shape[0])[:10] 	
#Plotting the partial dependence for the first feature in the dataset	
plot_pdp(rf_mine, X_train.iloc[idxs, :], y_train[idxs], X_train.columns[0])

Partial Dependence Plot for a particular feature

• 保持所有其他的列（特征的值）不变，在循环里面遍历考虑所有的单个特征，运行模型对于这个特征的每个独立的值预测结果。

• 绘制结果值，理解这个特征和目标变量背后的趋势。

保持其他的特征不变，改变依赖的特征，是解决外部噪声问题的关键，可以描述特征和目标的关系。

随机森林的树解释器

有一个非常有用的技术可以解释随机森林的流程，就是通过每个决策树的树解释器。

A very useful technique of explaining the flow of Random Forest predictions is through using Tree Interpreters for each of the Decision Trees.

class TreeInterpreter:	
    def predict(self, rf_model_tree, row):	
            prediction , bias, contribs = rf_model_tree.predict_row_for_ti(row, {})	
            print('prediction: '+str(prediction)+"bias: "+str(bias)+" contributions: "+str(contribs))	
            return prediction, bias, contribs

Tree Interpreter Class with the predict function

• 树的根节点的值（也就是所有训练数据的平均值）叫做偏差。

• 树解释器中的predict()函数调用了决策树中的predictrowfor_ti() 函数，预测了单个数据样本的结果。

• 对于每个分割，用来做分割的特征的分布使用分成子树之后样本的多样性的均值的数量来计算。

• 特征的分布存在一个数组中，直到叶子节点，结果返回给树解释器。

• 在叶子节点，样本的预测值返回，同时返回的还有偏差和贡献，也就是每个特征对于最终结果的贡献。

ti = TreeInterpreter()	
pr, bias, contribs  = ti.predict(rf_mine.trees[0], X_valid.iloc[0, :])	
#Waterfall Chart based on Tree Interpreter	
from waterfall import waterfall_chart	
waterfall_chart.plot(list(contribs.keys()), list(contribs.values()))

Function inside the Decision Tree, used by the Tree Interpreter

基于树解释器的瀑布图

每个特征到达叶子节点后的贡献可以通过瀑布图来可视化，使用python中的waterfall工具。

瀑布图是一个非常有效的方法来可视化每个特征的贡献，从根节点到每个叶子节点。

def predict_row_for_ti(self, row, feat_contribs):	
        if self.is_leaf: 	
           return self.value, self.parent_value , feat_contribs	
        if row[self.split_feature] < self.split_val:	
            if self.split_feature in feat_contribs.keys():	
                feat_contribs[self.split_feature] += self.leftTree.value - self.value	
            else:	
                feat_contribs[self.split_feature] = self.leftTree.value - self.value	
            return self.leftTree.predict_row_for_ti(row, feat_contribs)	
        else:	
            if self.split_feature in feat_contribs.keys():	
                feat_contribs[self.split_feature] += self.rightTree.value - self.value	
            else:	
                feat_contribs[self.split_feature] = self.rightTree.value - self.value	
            return self.rightTree.predict_row_for_ti(row, feat_contribs)

Waterfall Chart with the contributing features in X axis and resulting target values in the Y axis

功能增强

这个随机森林的库还有几个需要增强的地方：

• 构建树解释器，需要对随机森林做特征的相关性，找到频繁的进行分割的特征组，这些也许有高的相关性。

• 这个是第二种探索随机森林特征交互的方法，第一种方法是局部依赖绘图，可以扩展到超过一个特征。

注意：标准的机器学习库在速度和效率上是高度优化的，因此，你写的所有的代码需要通过Cython来运行，来获取更快的速度。你需要在代码的前面加上 ‘%%cython’ 。

最后的一些想法和本文重点

从头写一个机器学习的算法，这是一个很好的开始，可以帮助理解一些很有用的python库，包括Numpy和加速数学计算，Matplotlib进行数据可视化。

完整的实现在这里：https://github.com/SonaliDasgupta/MLandAIAlgorithmsFromScratch，代码进行了进一步的优化，我还在继续优化，欢迎大家提建议。

这对我是非常有趣的练习，我希望你可以对随机森林钻研的更深，理解这并不是什么复杂的“黑盒子”。

往期精彩回顾

1、最全的AI速查表|神经网络，机器学习，深度学习，大数据

2、深度学习论文阅读路线图

3、如何构建使用Python进行数据处理的肌肉记忆

4、Image-to-Image的论文汇总

5、我们从一阶段的物体检测器SSD，YOLOv3，FPN & Focal loss (RetinaNet）中学到了什么？

6、资源|10个机器学习和深度学习的必读免费课程

7、经验之谈|别再在CNN中使用Dropout了

8、我们从region based物体检测器 (Faster R-CNN, R-FCN, FPN)中能学到些什么？

9、非常好用的Python图像增强工具，适用多个框架

10、Kaggle竞赛介绍: Home Credit default risk（一）

本文可以任意转载，转载时请注明作者及原文地址。

请长按或扫描二维码关注我们