手撕机器学习算法 | 决策树 | DecisionTree

谭贤芳

已于 2024-07-16 18:30:25 修改

阅读量975

点赞数 23

分类专栏：手撕算法文章标签：机器学习算法决策树 c++ python numpy

于 2024-07-11 17:30:14 首次发布

本文链接：https://blog.csdn.net/m0_70254540/article/details/140281509

版权

手撕算法专栏收录该内容

7 篇文章 0 订阅

订阅专栏

❗写在前面 ❗：

… … … … … … … … . … … …
1️⃣ 提供两种语言实现：本文提供两种语言实现，分别为Python和C++，分别依赖的核心库为Numpy和Eigen.

2️⃣ 提供关键步的公式推导：写出详尽的公式推导太花时间，所以只提供了关键步骤的公式推导。更详尽的推导，可移步其他博文或相关书籍.

3️⃣ 重在逻辑实现：文中的代码主要以实现算法的逻辑为主，在细节优化上面，并没有做更多的考虑，如果发现有不完善的地方，欢迎指出.

🌞 欢迎关注专栏，相应的内容会持续更新，祝大家变得更强！！
… … … … … … … … . … … …

1、原理推导

简单来说，决策树可以看做是一组if-then规则的集合，为决策树的根结点到叶子结点的每一条路径都构建一条规则，路径中的内部结点特征代表规则条件，而叶子结点表示这条规则的结论。

完整的决策树模型包括特征选择、决策树构建和决策树剪枝三大部分。

为了构建一颗决策树，我们需要从训练集中不断选取具有良好分类能力的特征对数据集进行划分。那什么样的特征才算是“好特征”呢？

很简单，当通过某个特征划分之后，左右节点的类别尽可能的是同一类时，该特征就是一个“好特征”，因为它可以使划分后的左右节点“变纯”。

以上只是一个定性的解释，那如何定量得出哪个特征划分最优？这里可以用到三种常用的技术：信息增益、信息增益比、基尼系数。

信息熵

$E(D)=-\sum\limits_{k=1}^{Y}p_klogp_k$

其中 $D$ 表示数据集， $k$ 表示类别，且 $k$ ∈{1, 2, 3, … , $Y$ }， $p_k$ 表示该类别所占的比例。

熵是一种混乱程度，在这里可以理解为数据集中包含类别的“单纯”程度，当一个数据集所属类别越“单一”，其熵就越小，也可以说其混乱程度越低。

Python代码 —— 信息熵

import numpy as np

# 计算信息熵
def entropy(label: list):
    """
    label: 包含类别取值的列表
    返回计算得到的信息熵的值
    """
    # 类别去重
    uniquedLabel = set(label)
    
    # 统计每种类别的比例/概率
    probs = [label.count(lab) / len(label) for lab in uniquedLabel]

    # 计算信息熵
    entropyValue = - np.sum([pk * np.log2(pk) for pk in probs])

    return entropyValue

# 信息熵测试
data1 = [1, 2, 3, 4, 4, 4, 3, 2, 1]
data2 = [1, 1, 1, 2, 1, 2, 1, 2, 1]

print("data1的信息熵: ", entropy(data1))
print("data2的信息熵: ", entropy(data2))

# 输出：
data1的信息熵:  1.9749375012019272
data2的信息熵:  0.9182958340544896

可以直观的看出data2中的信息熵比data1中要小，我们的直观感受也是data2比data1更干净，更纯洁。

C++代码 —— 信息熵

#include <iostream>
#include <Eigen/Dense>
#include <algorithm>
#include <unordered_set>

// 计算信息熵
double entropy(const VectorXi& label)
{
    // 类别去重
    unordered_set<int> uniqueLabel(label.data(), label.data() + label.size());

    // 有Y种类别
    int Y = uniqueLabel.size();

    // 统计每种类别的概率并计算信息熵
    double entropyValue = 0.0;
    
    for (auto const& lab: uniqueLabel)
    {
        double p = static_cast<double>((label.array() == lab).count()) / label.size();
        entropyValue -= p * log2(p);
    }

    return entropyValue;
}

// 信息熵测试
int main()
{
    VectorXi data1(9);
    VectorXi data2(9);
    data1 << 1, 2, 3, 4, 4, 4, 3, 2, 1;
    data2 << 1, 1, 1, 2, 1, 2, 1, 2, 1;
    
    cout << "data1的信息熵为: " << entropy(data1) << endl;
    cout << "data2的信息熵为: " << entropy(data2) << endl;
    
    return 0;
}

// 输出：
data1的信息熵为: 1.97494
data2的信息熵为: 0.918296

信息增益

介绍了信息熵，信息增益就很好解释了。简单来说就是数据集 $D$ 在给定条件 $A$ 的情况下变为数据集 $D_A$ 所获得的不确定性增加：

$g (D, A) = E (D) - E (D ∣ A)$

用上面的写好的函数可以表示为：

$gain = entropy(D) - entropy(D_A)$

信息增益比

信息增益已经是一个很好的决定特征划分的指标了，但存在一个缺陷，就是当某个特征分类取值较多时，该特征的信息增益计算结果就会较大。所以，基于信息增益选择特征时，会偏向于取值较大的特征。而使用信息增益比可以解决这个问题：

$g_R(D, A)=\frac{g(D, A)}{E_A(D)}$

其中 $E_A(D)=-\sum\limits_{i=1}^{n}\frac{|D_i|}{D}log_2\frac{|D_i|}{D}$ ， $n$ 表示特征A的取值个数。

Python代码 —— 信息增益比

def gainRateDA(fi: ndarray, y: ndarray, A: float):
    """
	fi 某个特征数组
	y  标签数组
	A  划分的阈值条件
	"""
    # 这里的entropy函数就是上面实现的计算信息熵的
    ED = entropy(y.tolist())

    y1, y2 = y[fi < A], y[fi >= A]

    EDA =  len(y1) / len(y) * entropy(y1.tolist()) + \
           len(y2) / len(y) * entropy(y2.tolist())

    GDA = ED - EDA

    GRDA = GDA / (EDA + 1e-10)  # 添加一个非常小的数防止分母为0

    return GRDA

C++代码 —— 信息增益比

double gainRateDA(const VectorXd& fi, const VectorXi& y, double A)
{
	
    double ED = entropy(y);
    VectorXi y1(y.size()), y2(y.size());

    int n = 0, m = 0;
    for (int i = 0; i < fi.size(); i++)
    {
        if (fi(i) < A)
        {
            y1(n) = y(i);
            n++;
        }
        else
        {
            y2(m) = y(i);
            m++;
        }
    }
    y1.conservativeResize(n);
    y2.conservativeResize(m);
    double EDA = double(n) / y.size() * entropy(y1) + 
                 double(m) / y.size() * entropy(y2);
    double GDA = ED - EDA;
    double ERDA = GDA / (EDA + 1e-10);

    return ERDA;
}

基尼系数

基尼指数也是一种较好的特征选择方法。基尼指数是针对概率分布而言的。假设样本有 $K$ 个类，样本属于第 $k$ 类的概率为 $p_k$ ，则该样本类别概率分布的基尼指数可定义为：

$\sum\limits_{k=1}^{K}p_k(1-p_k)$

Python代码 —— 基尼指数

def gini(label):
   	"""
   	label 标签列表
   	"""
    probs = [label.count(lab) / len(label) for lab in set(label)]
    giniValue = np.sum([p * (1 - p) for p in probs])

    return giniValue

C++代码 —— 基尼指数

double gini(const VectorXi& label)
{	// 类别去重
    std::unordered_set<int> uniqueLabel(label.data(), label.data() + label.size());
   
    double giniValue = 0.0;
    for (auto const& lab: uniqueLabel)
    {
        double p = double((label.array() == lab).count()) / label.size();
        giniValue += p * (1 - p);
    }
    
    return giniValue;
}

以上三种特征选择方法，分别对应着三种经典的决策树算法：

{
	"ID3决策树": "信息增益",
	"C4.5决策树": "信息增益比",
	"CART决策树": "基尼指数"   # 既可以用来分类, 也可以用于回归, 下面也将重点对它进行实现
}

文字有点写不动了，了解到这里，我们直接上代码，用代码语言直接去感受决策树算法的原理，我尽量在代码中把注释写详细点。

2、算法实现

Python

# 这里先声明节点类型
# 该节点类型保存决策树构建时的各种信息以及子节点
import numpy as np
from numpy import ndarray, inf
from abc import ABC, abstractmethod
# 这里先声明节点类型
# 该节点类型保存决策树构建时的各种信息以及子节点
class TreeNode:
    def __init__(self, featrueIdx=None, threshold=None, leftBranch=None, rightBranch=None, leafValue=None):
        # 该节点划分的特征索引
        self.featrueIdx = featrueIdx
        # 该节点划分时的特征阈值
        self.threshold = threshold
        # 该节点的左支
        self.leftBranch = leftBranch
        # 该节点的右支
        self.rightBranch = rightBranch
        # 该节点的值
        self.leafValue = leafValue


# CART回归树和CART分类树, 基本上只有在计算不纯度和预测叶子节点的值方面不同
# 同时他们又都是基于二叉树结构的, 所以直接定义一个二叉决策树类BinaryDecisionTree作为基类
# 然后把calcImpurity, calcLeafValue两个方法作为抽象方法, 交给子类去实现
class BinaryDecisionTree(ABC):
    
    def __init__(self, minSampleNum=None, maxDepth=None, setMinImpurity=None):
        # 决策树的根节点
        self.root = None
        # 节点最少要包含的样本树
        self.minSampleNum = minSampleNum
        # 决策树构建时的最大深度
        self.maxDepth = maxDepth
        # 划分时最小不纯度值
        self.setMinImpurity = setMinImpurity
    
    # 留给子类去实现: 用于计算纯度的
    @abstractmethod
    def calcImpurity(self, y, y1, y2):
        pass

    # 留给子类去实现: 用于计算叶节点值的
    @abstractmethod
    def calcLeafValue(self, y):
        pass  
    
    # 划分数据集：根据阈值将数据集划分为左右两部分
    def _splitXy(self, Xy, fi, threshold):
        
        Xy_left = Xy[Xy[:, fi] < threshold]
        Xy_right = Xy[Xy[:, fi] >= threshold]

        return Xy_left, Xy_right

    # 递归构建决策树
    def constructTree(self, X, y, currentDepth=0):
        """
        X  训练的特征集
        y  对应的标签
        currentDepth  当前树的深度
        """
        
        # 设置当前不纯度为最大，在后续寻找划分特征和划分阈值时不断更新...
        minImpurity = inf
        # 最佳划分特征的索引
        bestFi = None
        # 最佳划分特征的阈值
        bestThreshold = None
        # 最佳划分条件下的左支X
        bestLeftBranch_X = None
        # 最佳划分条件下的左支y
        bestLeftBranch_y = None
        # 最佳划分条件下的右支X
        bestRightBranch_X = None
        # 最佳划分条件下的右支y
        bestRightBranch_y = None
        
        # 拼接合并一下Xy方便操作
        Xy = np.concatenate((X, y), axis=1)
        
        m, n = X.shape  # 样本数和特征维度

        # 如果当前样本数大于设置的最小样本数 且 当前深度未达最大深度时
        # 遍历寻找最佳划分特征及阈值, 条件是: 划分后的不纯度值最小 且
        # 还得小于 self.setMinImpurity
        if (m >= self.minSampleNum and currentDepth <= self.maxDepth):
            
            # 挨个取特征列
            for fi in range(n):
                
                fi_values = np.expand_dims(X[:, fi], axis=1)

                # 挨个取划分值
                for threshold in set(fi_values):

                    # 根据特征值划分数据集
                    Xy_left, Xy_right = self._splitXy(Xy, fi, threshold)
                    
                    # 划分后的数据集不能为空
                    if len(Xy_left != 0) and len(Xy_right != 0):
                        X_left = Xy_left[:, :n]
                        y_left = Xy_left[:, n:]
            
                        X_right = Xy_right[:, :n]
                        y_right = Xy_right[:, n:]

                        # 计算划分后的当前不纯度
                        currentImpurity = self.calcImpurity(y, y_left, y_right)
                        
                        # 判断当前的不纯度是否比最小不纯度低
                        # 若是就更新
                        if (currentImpurity < minImpurity):
                            bestFi = fi
                            bestLeftBranch_X = X_left
                            bestLeftBranch_y = y_left

                            bestRightBranch_X = X_right
                            bestRightBranch_y = y_right
                            bestThreshold = threshold
                            
                            minImpurity = currentImpurity
            
            # 经过上面两个for循环，我们找到了最佳划分参数
            # 如果当前最佳的划分参数中的不纯度的值，比我们设定的最小不纯度的值还要小
            # 则进行划分构造
            if (minImpurity < self.setMinImpurity):
                
                # 递归构造左支
                leftBranch = self.constructTree(bestLeftBranch_X, bestLeftBranch_y, currentDepth+1)
                # 递归构造右支
                rightBranch = self.constructTree(bestRightBranch_X, bestRightBranch_y, currentDepth+1)
                
                # 返回节点
                return TreeNode(featrueIdx=bestFi, 
                                threshold=bestThreshold, 
                                leftBranch=leftBranch, 
                                rightBranch=rightBranch, 
                                leafValue=None)
            
            # 不纯度已经满足要求就标记为叶节点
            # 并计算叶节点的值
            else:
                leafValue = self.calcLeafValue(y)
                return TreeNode(leafValue=leafValue)
            
    # 训练
    def fit(self, X, y):
        self.root = self.constructTree(X, y)

    # 预测单个样本的标签
    def _predict_x(self, x, treeNode=None):
        if treeNode is None:
            treeNode = self.root

        # 若是叶节点直接返回值
        if treeNode.leafValue is not None:
            return treeNode.leafValue
        
        # 判断当前值进入左节点还是右节点
        nextBranch = None
        featrueValue = x[treeNode.featrueIdx]
        if featrueValue < treeNode.threshold:
            nextBranch = treeNode.leftBranch
        else:
            nextBranch = treeNode.rightBranch
        
        # 走递归
        return self._predict_x(x, nextBranch)
    
    # 预测
    def predict(self, X):
        return [self._predict_x(x) for x in X]
    

def gini(label):
    
    probs = [label.count(lab) / len(label) for lab in set(label)]
    giniValue = np.sum([p * (1 - p) for p in probs])

    return giniValue


# CART分类树：实现一下BinaryDecisionTree中没有实现的两个方法即可
class ClassificationTree(BinaryDecisionTree):
    
    # 分类树计算不纯度的方法是算基尼系数
    def calcImpurity(self, y, y1, y2):
        p = len(y1) / len(y)
        
        return p * gini(y1) + (1 - p) * gini(y2)

    # 分类树决定叶节点的类别通过多数投票的方式
    def calcLeafValue(self, y):
        
        leafValue = None
        mode = 0
        for yi in set(y):
            yiNum= list(y).count(yi)
            if yiNum > mode:
                mode = yiNum
                leafValue = yi

        return leafValue        

# CART回归树：实现一下BinaryDecisionTree中没有实现的两个方法即可
class RegressionTree(BinaryDecisionTree):

    # 回归树计算不纯度的方法是方差减少量
    def calcImpurity(self, y, y1, y2):
        p = len(y1) / len(y)

        return np.sum(np.var(y, axis=0) - (p * np.var(y1, axis=0) + (1 - p) * np.var(y2, axis=0)))

    # 回归树计算节点值的方法是节点的均值
    def calcLeafValue(self, y):
        
        return np.mean(list(y))

下面的cpp代码就不写详细的注释了，只写关键步的逻辑注释，如果需要详细注释的话，在评论区说一下。

C++

#include <iostream>
#include <Eigen/Dense>
#include <vector>
#include <unordered_set>

using namespace Eigen;
using namespace std;

// 树节点
struct TreeNode
{
    double threshold;
    int featrueIdx;
    TreeNode* leftBranch;
    TreeNode* rightBranch;
    double leafValue;

    TreeNode() : threshold(0.0), featrueIdx(0), leftBranch(nullptr), rightBranch(nullptr), leafValue(0.0) {}
};

// 二叉决策树基类
class BinaryDecisionTree
{
private:
    TreeNode* m_root;
    int m_minSampleNum;
    int m_maxDepth;
    double m_setMinImpurity;

public:
    BinaryDecisionTree(int minSampleNum, int maxDepth, double setMinImpurity)
        : m_root(nullptr), m_minSampleNum(minSampleNum), m_maxDepth(maxDepth), m_setMinImpurity(setMinImpurity)
    {
    }

    // 计算不纯度
    virtual double calcImpurity(const VectorXd& y, const VectorXd& y1, const VectorXd& y2) = 0;

    // 计算叶子值
    virtual double calcLeafvalue(const VectorXd& y) = 0;

    // 构造决策树
    TreeNode* construct(const MatrixXd& X, const VectorXd& y, int currentDepth = 0)
    {
        double minImpurity = std::numeric_limits<double>::infinity();
        int bestFi = 0;
        double bestThreshold = 0.0;

        MatrixXd bestLeftBranch_X;
        VectorXd bestLeftBranch_y;
        MatrixXd bestRightBranch_X;
        VectorXd bestRightBranch_y;

        int m = X.rows();
        int n = X.cols();

        if (m > m_minSampleNum && currentDepth < m_maxDepth)
        {
            for (int fi = 0; fi < n; fi++)
            {
                VectorXd fi_values = X.col(fi);
                std::unordered_set<double> fi_values_set(fi_values.data(), fi_values.data() + fi_values.size());

                for (const double& threshold : fi_values_set)
                {
                    std::vector<int> a, b;
                    for (int idx = 0; idx < m; idx++)
                    {
                        if (fi_values[idx] < threshold)
                            a.push_back(idx);
                        else
                            b.push_back(idx);
                    }

                    if (a.size() != 0 && b.size() != 0)
                    {
                        MatrixXd X_left = X(Eigen::Index(a.size()), Eigen::all);
                        MatrixXd X_right = X(Eigen::Index(b.size()), Eigen::all);
                        VectorXd y_left = y(Eigen::Index(a.size()));
                        VectorXd y_right = y(Eigen::Index(b.size()));

                        double currentImpurity = this->calcImpurity(y, y_left, y_right);

                        if (currentImpurity < minImpurity)
                        {
                            bestFi = fi;
                            bestLeftBranch_X = X_left;
                            bestLeftBranch_y = y_left;
                            bestRightBranch_X = X_right;
                            bestRightBranch_y = y_right;
                            bestThreshold = threshold;

                            minImpurity = currentImpurity;
                        }
                    }
                }
            }

            if (minImpurity < m_setMinImpurity)
            {
                TreeNode* leftBranch = construct(bestLeftBranch_X, bestLeftBranch_y, currentDepth + 1);
                TreeNode* rightBranch = construct(bestRightBranch_X, bestRightBranch_y, currentDepth + 1);

                TreeNode* Node = new TreeNode();
                Node->featrueIdx = bestFi;
                Node->threshold = bestThreshold;
                Node->leftBranch = leftBranch;
                Node->rightBranch = rightBranch;

                return Node;
            }
            else
            {
                double leafValue = calcLeafvalue(y);
                TreeNode* Node = new TreeNode();
                Node->leafValue = leafValue;

                return Node;
            }
        }

        return nullptr;
    }

    // 训练
    void fit(const MatrixXd& X, const VectorXd& y)
    {
        if (m_root != nullptr)
        {
            deleteTree(m_root); // 删除当前树
        }
        m_root = construct(X, y);
    }

    // 删除树的函数，辅助函数
    void deleteTree(TreeNode* node)
    {
        if (node == nullptr)
            return;
        deleteTree(node->leftBranch);
        deleteTree(node->rightBranch);
        delete node;
    }

    // 预测单样本
    double predict_x(const RowVectorXd& x, const TreeNode* treeNode = nullptr)
    {
        if (treeNode == nullptr)
            treeNode = m_root;

        if (treeNode->leafValue != 0.0)
            return treeNode->leafValue;

        TreeNode* nextNode = nullptr;
        double featureValue = x(treeNode->featrueIdx);

        if (featureValue < treeNode->threshold)
            nextNode = treeNode->leftBranch;
        else
            nextNode = treeNode->rightBranch;

        return predict_x(x, nextNode);
    }

    // 预测全部样本
    VectorXd predict(const MatrixXd& X)
    {
        VectorXd y_pred(X.rows());

        for (int i = 0; i < X.rows(); ++i)
        {
            y_pred(i) = predict_x(X.row(i));
        }

        return y_pred;
    }
};
// 计算基尼不纯度
double gini(const VectorXd& label)
{
    std::unordered_set<int> uniqueLabel(label.data(), label.data() + label.size());
    
    double giniValue = 0.0;
    for (auto const& lab : uniqueLabel)
    {
        double p = double((label.array() == lab).count()) / label.size();
        giniValue += p * (1 - p);
    }
    
    return giniValue;
}

// 分类决策树类
class ClassificationTree : public BinaryDecisionTree
{
public:
    ClassificationTree(int minSampleNum, int maxDepth, double setMinImpurity)
        : BinaryDecisionTree(minSampleNum, maxDepth, setMinImpurity)
    {
    }

    // 计算分类决策树的不纯度
    double calcImpurity(const VectorXd& y, const VectorXd& y1, const VectorXd& y2) override
    {
        double p = double(y1.size()) / y.size();
        return p * gini(y1) + (1 - p) * gini(y2);
    }

    // 计算叶子值
    double calcLeafvalue(const VectorXd& y) override
    {
        double leafValue = 0.0;
        int mode = 0;
        std::unordered_set<double> unique_y(y.data(), y.data() + y.size());

        for (const auto& yi : unique_y)
        {
            int yiNum = (y.array() == yi).count();
            if (yiNum > mode)
            {
                mode = yiNum;
                leafValue = yi;
            }
        }
        return leafValue;
    }
};

// 计算方差
double variance(const VectorXd& y)
{
    return (y.array() - y.mean()).array().square().sum() / y.size();
}

// 回归决策树类
class RegressionTree : public BinaryDecisionTree
{
public:
    RegressionTree(int minSampleNum, int maxDepth, double setMinImpurity)
        : BinaryDecisionTree(minSampleNum, maxDepth, setMinImpurity)
    {
    }

    // 计算回归决策树的不纯度
    double calcImpurity(const VectorXd& y, const VectorXd& y1, const VectorXd& y2) override
    {
        double p = double(y1.size()) / y.size();
        return variance(y) - (p * variance(y1) + (1 - p) * variance(y2));
    }

    // 计算叶子值
    double calcLeafvalue(const VectorXd& y) override
    {
        return y.mean();
    }
};