python程序实现CART回归树算法

算法原理

CART树属于决策树模型，它可以用于分类和回归问题，这两者的不同主要体现在特征选择上，CART分类树基于Gini系数最小化，而CART回归树是基于误差平方和最小化。另外，CART树区别于其他树模型的特点之一是它是一类二叉树模型，每次进行特征选择时都将数据集分成“是”或“否”这两种。本文将讲解CART回归树模型的原理以及程序实现。

模型

已知一组包含$n$个样本的训练集
{ x i , y i } i = 1 n , x i ∈ R t , y i ∈ R \{x_i,y_i\}_ {i=1}^n,x_i\in\mathbb R^t,y_i\in\mathbb R {xi​,yi​}i=1n​,xi​∈Rt,yi​∈R
其中，每个样本$x_i$都有$t$个特征（假设每个特征都是连续值），输出变量$y_i$也是连续值。现在的目的是希望通过对训练集进行训练，生成决策树。之后往模型中代入新的样本$x_j$，得到预测的连续值$y_j$。

原理

CART树生成主要包括
1、遍历选择最优特征
2、按特征进行分类
3、迭代前两步
4、停止迭代

遍历选择最优特征

CART回归树的特征选择是基于误差平方和最小化。以一个节点的最优特征选择为例，因为知道一个节点是如何划分的也就可以知道如何划分所有节点，这是一个迭代的过程。假设当前节点的数据集为$N$，依次将每个特征的每个取值$s$作为阈值，将数据集分为两类，分别为
N 1 = { ( x , y ) ∣ A j ( x ) ≤ s , ( x , y ) ∈ N } N_1=\{(x,y)|A_j(x)\leq{s},(x,y)\in{N}\} N1​={(x,y)∣Aj​(x)≤s,(x,y)∈N}
N 2 = { ( x , y ) ∣ A j ( x ) > s , ( x , y ) ∈ N } N_2=\{(x,y)|A_j(x)>{s},(x,y)\in{N}\} N2​={(x,y)∣Aj​(x)>s,(x,y)∈N}
其中，$A_j(x)$表示样本$x$的第$j$个特征的取值，得到这两类数据集后，便可以计算它们的预测输出值$\hat{y}$，是它们各自的数据集中$y$变量的均值。它们的预测输出值分别为
$${\hat{y}}_{N_1}=\frac{1}{|N_1|}\sum _{(x,y)\in N_1}y$$
$${\hat{y}}_{N_2}=\frac{1}{|N_2|}\sum _{(x,y)\in N_2}y$$
然后计算这两类的误差平方总和
$$e(N,j,s)=\sum _{(x,y)\in N_1}(y-{\hat{y}}_{N_1}{)}^2+\sum _{(x,y)\in N_2}(y-{\hat{y}}_{N_2}{)}^2$$
因为数据集$N$的样本是有限的，样本的特征是有限的，特征的取值也是有限的，因此可以通过遍历所有特征的所有取值，得到一个$e(N,j,s)$序列，并从该序列中选择数值最小的一个作为当前节点的特征划分。

得到最优的特征及对应的取值之后便可以将数据集分成两个小数据集，也即是两个子节点，然后用与上面同样的操作继续对子节点进行分类。

停止迭代

连续的特征与离散特征在停止迭代的条件上有所不同。

对于离散的特征来说，对某个特征进行分类，分类之后得到的子数据组中该特征的取值都是相同的。因此在后续的分类中，不必再选取该特征进行分类。所以随着分类次数的增多，可选的分类特征会越来越少，当没有可选的特征时便可停止迭代。以西瓜数据集为例，当选择到分类特征为"色泽"时，节点数据集将分成三个子数据集，分别是"青绿"，“乌黑"和"浅白”，每个子数据集中的数据样本在特征"色泽"上的取值都是一样的，因此之后的分类不需要再使用该特征。

而对于连续的特征来说，每次都是以某个特征的某个取值作为阈值进行二分类，这样分类之后得到的小数据组仍可以继续对该特征进行分类。所以此时的迭代停止条件与上面不同，可以是节点样本数低于某个设定的阈值或树的深度高于某个设定的阈值。

程序实现

各函数

CART回归树的生成包含3个函数，分别是计算误差平方和，选择最优特征及取值，树生成主函数。

计算误差平方和

def cal_mse(y):
    """误差平方和计算

    参数
    -------
    y:真实值, 类型:narray, shape:{n_samples}

    返回
    ----
    m:误差平方和, 类型:float
    """
    return y.var()*y.size

在训练时，某个节点的误差平方和是指实际值与预测值的残差平方和，处于同一个节点上的样本的预测都是相同的，假设实际值为$y_i$，$(i=1,2,...,s)$，则预测值为$\frac{1}{s}{\sum }_{i=1}^s{y}_i$，那么误差平方和为
$$\sum _{j=1}^s(y_j-\frac{1}{s}\sum _{i=1}^s{y}_i{)}^2=s\ast VAR(y_j)$$

选择最优特征及取值

def choose_features(X, y):
    """选择连续特征及对应的取值（单个）

    参数
    ------
    X:特征, 类型:ndarray, shape:{n_samples, n_features}

    y:类别, 类型:ndarray, shape:{n_samples}

    返回
    -----
    fea_index:选出的特征, 类型:integer

    fea_c:特征阈值, 类型:float

    min_mse:最小误差平方和, 类型:float
    """
    n_samples, n_features = X.shape

    # 最优特征的索引
    fea_index = 0
    # 最优特征的阈值
    fea_c = 0
    # 最小误差平方和
    min_mse = float('inf')

    for i in range(n_features):
        # 对第i个特征的数值从大到小进行排序
        X_value = sorted(set(X[:,i]))
        for j in range(len(X_value)-1):
            cur_c = np.mean(X_value[j:j+2]) # 当前阈值，在当前索引j与j+1之间进行划分，小于等于阈值的一组，大于阈值的另一组
            cur_mse = cal_mse(y[X[:,i]<=cur_c])+cal_mse(y[X[:,i]>cur_c])
            if cur_mse<min_mse:
                min_mse = cur_mse
                fea_index = i
                fea_c = cur_c
    return fea_index, fea_c, min_mse

构建CART回归树

def tree_regress(X, y, samples=10):
    """构建决策树

    参数
    ------
    X:特征, 类型:ndarray, shape:{n_samples, n_features}

    y:类别编号, 类型:ndarray, shape:{n_samples}

    samples:样本数限制, 类型:integer
    
    返回
    -----
    dic:决策树, 类型:dict
    """
    if not len(X):return 
    # 样本数限制
    if len(X)<=samples:return np.mean(y)

    n_samples, n_features = X.shape

    fea_index,fea_c = choose_features(X, y)[:2]
    # 决策树构建 key：left表示小于等于特征值，right表示大于特征值
    dic = {}
    dic['fea_name'] = X_name[fea_index] # 特征名
    dic['fea_c'] = fea_c # 特征取值
    fea_value = X[:, fea_index]
    bool_l = fea_value<=fea_c
    bool_r = fea_value>fea_c
    dic['left'] = tree_regress(X[bool_l], y[bool_l])
    dic['right'] = tree_regress(X[bool_r], y[bool_r])
    return dic 

实例化演示

# 导库
import numpy as np
from sklearn.datasets import load_boston
# 训练数据
data = load_boston()
X = data.data
y = data.target
X_name = data.feature_names
# 生成树
tree_regress(X, y)

因为没有进行后剪枝，所以得到的树结构过于累赘，这里不做展示。

封装成一个类

增加了训练函数fit和预测函数predict

class Tree_Regress:
    def __init__(self, samples=10):
        self.samples = samples # 最小样本数

    def cal_mse(self, y):
        """误差平方和计算

        参数
        -------
        y:真实值, 类型:narray, shape:{n_samples}

        返回
        ----
        m:误差平方和, 类型:float
        """
        return y.var()*y.size
    
    def choose_features(self, X, y):
        """选择连续特征及对应的取值（单个）

        参数
        ------
        X:特征, 类型:ndarray, shape:{n_samples, n_features}

        y:类别, 类型:ndarray, shape:{n_samples}

        返回
        -----
        fea_index:选出的特征, 类型:integer

        fea_c:特征阈值, 类型:float

        min_mse:最小误差平方和, 类型:float
        """
        n_samples, n_features = X.shape

        # 最优特征的索引
        fea_index = 0
        # 最优特征的阈值
        fea_c = 0
        # 最小误差平方和
        min_mse = float('inf')

        for i in range(n_features):
            # 对第i个特征的数值从大到小进行排序
            X_value = sorted(set(X[:,i]))
            for j in range(len(X_value)-1):
                cur_c = np.mean(X_value[j:j+2]) # 当前阈值，在当前索引j与j+1之间进行划分，小于等于阈值的一组，大于阈值的另一组
                cur_mse = self.cal_mse(y[X[:,i]<=cur_c])+self.cal_mse(y[X[:,i]>cur_c])
                if cur_mse<min_mse:
                    min_mse = cur_mse
                    fea_index = i
                    fea_c = cur_c
        return fea_index, fea_c, min_mse

    def tree_regress(self, X, y):
        """构建决策树

        参数
        ------
        X:特征, 类型:ndarray, shape:{n_samples, n_features}

        y:类别编号, 类型:ndarray, shape:{n_samples}

        返回
        -----
        dic:决策树, 类型:dict
        """
        if not len(X):return 
        # 样本数限制
        if len(X)<=self.samples:return np.mean(y)

        n_samples, n_features = X.shape
        
        fea_index,fea_c = self.choose_features(X, y)[:2]
        # 决策树构建 key：left表示小于等于特征值，right表示大于特征值
        dic = {}
        dic['fea_name'] = X_name[fea_index]
        dic['fea_c'] = fea_c
        fea_value = X[:, fea_index]
        bool_l = fea_value<=fea_c
        bool_r = fea_value>fea_c
        dic['left'] = self.tree_regress(X[bool_l], y[bool_l])
        dic['right'] = self.tree_regress(X[bool_r], y[bool_r])
        return dic 

    def check(self, tree, X):
        """预测
        """
        if not tree or not len(X):return
        if isinstance(tree, float):return tree # 判断是否遍历到了叶子节点
        cur_fea_name = tree['fea_name']
        cur_fea_index = np.where(self.X_name==cur_fea_name)[0][0]
        cur_fea_c = tree['fea_c']
        
        if X[cur_fea_index]<=cur_fea_c:
            return self.check(tree['left'], X)
        else:
            return self.check(tree['right'], X)

    def fit(self, X, y, X_name):
        self.X_name = X_name
        self.tree = self.tree_regress(X, y)

    def predict(self, X):
        res = []
        for i in range(len(X)):
            res.append(self.check(self.tree, X[i]))
        return np.array(res)

演示

clf = Tree_Regress()
clf.fit(X, y, X_name)
predict_y = clf.predict(X)

不足

没有对生成的树进行必要的后剪枝，导致模型结构过于复杂，增加了过拟合的风险，关于剪枝的操作将在以后补充。

----end----

参考资料
李航《统计学习方法》