决策树模型的各种Ensemble

概览

Ensemble的方法主要有两大类:BaggingBoosting

  1. Boosting主要关注降低偏差,因此Boost能基于泛化性能相当弱的学习器构建出很强的集成;
  2. Bagging主要关注降低方差,因此它在不剪枝的决策树、神经网络等学习器上效用更为明显。
  3. Boosting的个体学习器之间存在强依赖关系,必须串行生成;
  4. Bagging的个体学习器之间不存在强依赖关系,可以同时生成即并行化

Bagging

基础Bagging

先讲Bagging,Bagging是Bootstrap aggregation的缩写。所谓的Bootstrap有放回抽样,而这里的抽样指的是对数据样本的抽样。

如果对于一个有$n$个样本的数据集$D$,有放回抽样$n’$个数据,那么当$n$足够大的时候,满足抽样出来的样本个数(无重复的个数)和原数据集的比例是:$(1 - 1/e) (≈63.2\%)$

证明

某个样本没有被抽中的概率是:$p_{not} = (1-\frac{1}{n})^n$
$\frac{1}{p_{not}} = (\frac{n}{n-1})^{n} = (1+\frac{1}{n-1})^{n-1}(1+\frac{1}{n-1})$

当n很大时,上式等于e(根据常用极限:$lim_{x\rightarrow∞}(1+\frac{1}{x})^x=e$)。

因此,$p_{not} = (1-\frac{1}{n})^n = \frac{1}{e}$。
回到Bagging,Bagging的基本做法

  • 1 从样本有放回地抽取n个样本;
  • 2 在所有的属性上,对这n个样本建立分类器;
  • 3 重复上述过程m次,得到m个分类器;
  • 4 将数据放在这m个分类器上分类,最终结果由所有分类器结果投票决定。

Bagging实践

sklearn.ensemble.BaggingClassifier提供了Bagging的模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np  
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.ensemble import BaggingClassifier
from sklearn import datasets

#读取数据,划分训练集和测试集
iris=datasets.load_iris()
x=iris.data[: , :2]
y=iris.target
x_train, x_test, y_train, y_test = train_test_split(x, y, train_size=0.7, random_state=1)

#模型训练
# sklearn 自带的决策树分类器
model1=DecisionTreeClassifier(max_depth=3)
# sklearn自带的bagging分类器
model2=BaggingClassifier(model1,n_estimators=100,max_samples=0.3)
model1.fit(x_train,y_train)
model2.fit(x_train,y_train)
model1_pre=model1.predict(x_test)
model2_pre=model2.predict(x_test)
res1=model1_pre==y_test
res2=model2_pre==y_test
print '决策树测试集正确率%.2f%%'%np.mean(res1*100)
print 'Bagging测试集正确率%.2f%%'%np.mean(res2*100)
1
2
决策树测试集正确率75.56%
Bagging测试集正确率77.78%

随机森林

随机森林(Random Forest, 简称RF)是Bagging的一个扩展变体。

RF在Bagging的基础上,加入了随机属性选择,即,对特征进行无放回抽样,并且使用CART树。

随机森林的基本做法

  • 1 首先在样本集中有放回的抽样n个样本;
  • 2 在所有的属性当中再随机选择K个属性;
  • 3 根据这n个样本的K个属性,建立CART树;
  • 4 重复以上过程m次,得到了m棵CART树;
  • 5 利用这m棵树对样本进行预测并投票。

随机森林实践

sklearn.ensemble.RandomForestClassifier提供了随机森林的模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np  
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn import datasets

#读取数据,划分训练集和测试集
iris=datasets.load_iris()
x=iris.data[:,:2]
y=iris.target
x_train, x_test, y_train, y_test = train_test_split(x, y, train_size=0.7, random_state=1)

#模型训练
model1=DecisionTreeClassifier(max_depth=3)
model2=RandomForestClassifier(n_estimators=200, criterion='entropy', max_depth=3)
model1.fit(x_train,y_train)
model2.fit(x_train,y_train)
model1_pre=model1.predict(x_train)
model2_pre=model2.predict(x_train)
res1=model1_pre==y_train
res2=model2_pre==y_train
print '决策树训练集正确率%.2f%%'%np.mean(res1*100)
print '随机森林训练集正确率%.2f%%'%np.mean(res2*100)
1
2
决策树训练集正确率83.81%
随机森林训练集正确率85.71%

Boosting

Boosting(提升)通过给样本设置不同的权值,每轮迭代调整权值。

不同的提升算法之间的差别,一般是:

  1. 如何更新样本的权值
  2. 如何组合每个分类器的预测,即,调整分类器的权值
    其中Adaboost中,样本权值是增加那些被错误分类的样本的权值,分类器$C_i$的重要性依赖于它的错误率。

AdaBoost

直接看算法:

  • 1.初始化训练数据的权值分布:$D_1 = (w_{11}, …, w_{1i}, …, w_{1N},)$,$w_{1,i}=\frac{1}{N}$。(N是数据个数)
  • 对$m=1, …, M$(M是弱分类器的个数):
    • 2.使用具有权值分布$D_m$训练得到弱分类器:$G_m(x)$
    • 3.计算$G_m(x)$在训练集上的分类错误率:$e_m =\sum_{i=1}^Nw_{mi}I(G_m(x_i)≠y_i)$
    • 4.计算$G_m(x)$的系数:$α_m=\frac{1}{2}log\frac{1-e_m}{e_m}$(分类器错误率越大,权重越小
    • 5.更新训练集权重:$w_{m+1, i} = \frac{w_{mi}}{Z_m}exp(-α_my_iG_m(x_i))$(**$x_i$分类错误则提高权重**)
  • 6.构建弱分类器的线性组合:$f(x) = \sum_{m=1}^Mα_mG_m(x)$

GBDT

GBDT(Gradient Boosting Decision Tree),梯度提升算法又叫 **MART(Multiple Additive Regression Tree)**。

从后面的名字可以看出,这是一种回归树的模型。它使用的基础分类器也是CART。

着重分析GBDT中的Boosting,即Additive Training

这里的方法是:第一轮去拟合一个大致上和目标差不多的值,然后计算残差,下一轮的拟合目标就是这个残差,即:

  • 1.用训练集训练一个弱分类器:$f_1(x) = y$
  • 2.用残差训练一个弱分类器:$h_1(x) = y - f_1(x)$
  • 3.获得新模型:$f_2(x) = f_1(x) + h_1(x)$ 上一轮迭代得到的强学习器是$f_{t-1}(x)$,损失函数是$L(y, f_{t-1}(x))$。

我们本轮迭代的目标是:找到一个CART回归树模型的弱学习器$h_t(x)$,让本轮的损失损失$L(y, f_{t}(x)) =L(y, f_{t-1}(x)+ h_t(x))$最小。

  • 第t轮的第i个样本,损失函数的负梯度表示:$r_{ti} = -\bigg[\frac{\partial L(y,f(x_i)))}{\partial f(x_i)}\bigg]_{f(x) = f_{t-1}(x)}$,负梯度可以用来估计残差(注:为什么要使用负梯度而不直接使用残差是为了使用不同Loss function时,使用负梯度更容易优化)。
  • 再用$(x_i, r_{ti})$去拟合下一棵CART回归树:
    • 对$j=1,2,…,J$,计算:$c_{tj} = \underbrace{arg; min}_{c}\sum\limits_{x_i \in R_{tj}} L(y_i,f_{t-1}(x_i) +c)$
  • 得到本轮的拟合函数:$h_t(x) = \sum\limits_{j=1}^{J}c_{tj}I(x \in R_{tj})$
    • 更新强学习器:$f_{t}(x) = f_{t-1}(x) + \sum\limits_{j=1}^{J}c_{tj}I(x \in R_{tj})$

上面算法描述版本参考自李航老师的《统计学习方法》,$c_{tj}$是决策树单元$R_{tj}$上固定输出值。然而更常见的算法版本(参考wikipedia)并不使用$c$,而是使用**步长$\gamma$**,步长×负梯度更reasonable:

即:

  • 使用$(x_i, r_{ti})$去拟合下一棵CART回归树:
    • 拟合的输出是$h_m(x)$
      -再用线性搜索找到最佳步长$\gamma_m = \underset{\gamma}{arg; min}\sum_{i=1}^nL(y_i, F_{m-1}(x_i)+\gamma h_m(x_i))$,其中的$\gamma h_m(x)$等价于$c$。

XGBoost

XGBoost的符号表示稍微有些不同,参考陈天奇slide

XGBoost本质上也是一个gradient boosting。

我们每轮都训练一个基础分类器:

$\hat y_i^{(0)} = 0$

$\hat y_i^{(1)} = f_1(x_i) = \hat y_i^{(0)} + f_1(x_i)$

$\hat y_i^{(2)} = f_1(x_i) + f_2(x_i) = \hat y_i^{(1)} + f_2(x_i)$

$\hat y_i^{(t)} = \sum_{k=1}^t f_k(x_i) = \hat y_i^{(t-1)} + f_t(x_i)$

上面的最后一条式子解释如下:当我们训练了t轮以后,有了t个弱学习器,每一轮都将之前已经Boosting的强学习器和现在这轮的弱学习器的结果相加。

目标函数(损失函数):$Obj^{(t)} = \sum_{i=1}^nl(y_i, \hat y_i^{(t)})+\sum_{i=1}^tΩ(f_i)$

回忆泰勒展开:$f(x + \Delta x)≈f(x)+f’(x)\Delta x+\frac{1}{2}f’’(x)\Delta x^2$

将$\hat y_i^{(t)} = \hat y_i^{(t-1)} + f_t(x_i)$带入目标函数,转换成:

$Obj^{(t)} = \sum_{i=1}^n[l(y_i, \hat y_i^{(t-1)})+g_if_t(x_i)+\frac{1}{2}h_if_t^2(x_i)]+Ω(f_t)+constant$,其中,$g_i = \partial_{\hat y_i^{(t-1)}}l(y_i, \hat y_i^{(t-1)})$,$h_i = \partial^2_{\hat y_i^{(t-1)}}l(y_i, \hat y_i^{(t-1)})^2$

把上面的常数都去掉,剩下:

$Obj^{(t)} = \sum_{i=1}^n[g_if_t(x_i)+\frac{1}{2}h_if_t^2(x_i)]+Ω(f_t)$

引入叶结点的权重:$f_t(x) = w_{q(x)}$,$q(x)$是样本到叶结点的映射,正则函数:$Ω(f_t) = \gamma T + \frac{1}{2}\lambda \sum^T_{j=1}w^2_j$。$T$是叶结点的个数。

定义在叶结点j的样本集:$I_j = {i|q(x_i) = j}$
把叶结点权重和正则函数带入目标函数,得到:

$$Obj^{(t)} = \sum_{i=1}^n[g_iw_{q(x_i)}+\frac{1}{2}h_iw^2_{q(x_i)}]+\gamma T +
\frac{1}{2}\lambda \sum^T_{j=1}w^2_j$$

$$=\sum_{j=1}^T[(\sum_{i∈I_j }g_i)w_j +
\frac{1}{2}(\sum_{i∈I_j}h_i+\lambda)w_j^2]+\gamma T$$

定义:$G_j = \sum_{i∈I_j}g_i$,$H_j = \sum_{i∈I_j}h_i$,

$Obj^{(t)} =\sum_{j=1}^T[G_jw_j + \frac{1}{2}(H_j + \lambda)w_j^2] + \lambda T$

回忆一下一元二次函数的性质:

对于:$Gx + \frac{1}{2}Hx^2$,($H>0$),最小值为:$-\frac{1}{2}\frac{G^2}{H}$,在$ -\frac{G}{H}$处取得。

回到目标函数中去,如果树的结构($q(x)$)固定,那最优的权值分配是:

$w_j^* = -\frac{G_i}{H_j+\lambda}$

$Obj^* = -\frac{1}{2}\sum^T_{j=1}\frac{G_j^2}{H_j+\lambda}+\gamma T$

接下来考虑如何学习树的结构(Greedy Learning)

回忆一下CART是如何做的:

对现有特征A的每一个特征,每一个可能的取值a,根据样本点对$A=a$的测试是“是”还是“否”,将$D$分割成$D_1$和$D_2$两部分,计算$A=a$时的基尼指数。选择基尼指数最小的特征机器对应的切分点作为最优特征最优切分点

这里不算基尼指数,这里的Gain是:$\frac{1}{2}[\frac{G^2_L}{H_L+\lambda}+\frac{G^2_R}{H_R+\lambda}-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}] - \gamma$。

上式方括号中的三项分别代表:左子树的得分右子树的得分如果不分割的得分。得分可以理解成是损失的反面。
寻找最优分割的算法:

  • 对每个结点,遍历所有feature:
    • 对每个feature,将数据样本按照feature 值排序;
    • 使用linear scan决定这个feature的最佳split;
    • 如此循环,找到所有feature的最佳split。

上面算法的时间复杂度是:$O(d\cdot K\cdot nlogn )$,其中,$n$是数据个数,$d$是特征个数,$K$是树深度。

剪枝

这里的剪枝和其他普通的决策树剪枝没有区别,分为前剪枝和后剪枝。其中,后剪枝的策略是:递归地减掉所有产生negative gain的split。

Boosting实践

sklearn中自带:sklearn.ensemble.GradientBoostingClassifiersklearn.ensemble.AdaBoostClassifier

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np  
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import AdaBoostClassifier
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn import datasets

iris=datasets.load_iris()
x=iris.data[:,:2]
y=iris.target

model1=DecisionTreeClassifier(max_depth=5)
model2=GradientBoostingClassifier(n_estimators=100)
model3=AdaBoostClassifier(model1,n_estimators=100)
model1.fit(x,y)
model2.fit(x,y)
model3.fit(x,y)
model1_pre=model1.predict(x)
model2_pre=model2.predict(x)
model3_pre=model3.predict(x)
res1=model1_pre==y
res2=model2_pre==y
res3=model3_pre==y
print '决策树正确率%.2f%%'%np.mean(res1*100)
print 'GDBT正确率%.2f%%'%np.mean(res2*100)
print 'AdaBoost正确率%.2f%%'%np.mean(res3*100)
1
2
3
决策树正确率84.67%
GDBT正确率92.00%
AdaBoost正确率92.67%

xgboost实践

xgboost需要额外安装一下:官方安装地址

这里简单说下ubuntu下python接口的安装:

1
2
3
4
5
git clone --recursive
https://github.com/dmlc/xgboost
cd xgboost; make -j4
cd python-package; sudo
python setup.py install

同样,官网这里对调参方法有很详细的介绍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import xgboost as xgb  
from sklearn.model_selection import train_test_split
from sklearn import datasets

iris=datasets.load_iris()
x=iris.data[:,:2]
y=iris.target
x_train, x_test, y_train, y_test = train_test_split(x, y, train_size=0.7, random_state=1)
data_train = xgb.DMatrix(x_train,label=y_train)
data_test=xgb.DMatrix(x_test,label=y_test)
param = {}
param['objective'] = 'multi:softmax'
param['eta'] = 0.1
param['max_depth'] = 5
param['silent'] = 1
param['nthread'] = 4
param['num_class'] = 3
watchlist = [ (data_train,'train'), (data_test, 'test') ]
num_round = 10
bst = xgb.train(param, data_train, num_round, watchlist );
pred = bst.predict( data_test );
print ('predicting, classification error=%f' % (sum( int(pred[i]) != y_test[i] for i in range(len(y_test))) / float(len(y_test)) ))
1
2
3
4
5
6
7
8
9
10
11
[0]	train-merror:0.133333	test-merror:0.266667
[1] train-merror:0.142857 test-merror:0.266667
[2] train-merror:0.12381 test-merror:0.266667
[3] train-merror:0.12381 test-merror:0.266667
[4] train-merror:0.12381 test-merror:0.266667
[5] train-merror:0.114286 test-merror:0.288889
[6] train-merror:0.104762 test-merror:0.288889
[7] train-merror:0.104762 test-merror:0.288889
[8] train-merror:0.114286 test-merror:0.288889
[9] train-merror:0.114286 test-merror:0.288889
predicting, classification error=0.288889

参考链接