Machine Learning-决策树与集成学习

决策树与集成学习

目录

·简介

·决策树

·booststrap

·bagging

·boosting

·随机森林

·Adaboost

·GBDT

·XGboost

简介

本文主要针对集成学习中的随机森林、Adaboost和GBDT做详细介绍,因为介绍这两个模型,必须要从最基础的决策树讲起,所以对决策树和一些集成方法也做了详细介绍。所以本文的框架是从决策树讲起,再过渡到集成方法,再过度到最好的随机森林、Adaboost和GBDT。

决策树

决策树(Decision Tree)是一种非参数的有监督学习方法,它能够从一系列有特征和标签的数据中总结出决策规则,并用树状图的结构来呈现这些规则,以解决分类和回归问题。
优点:简单易懂,容易解释,可视化,适用性广
缺点:容易过拟合,数据中的小变化会影响结果,不稳定,每一个节点的选择都是贪心算法,不能保证全局最优解。

决策树的组成:

根节点:第一个选择点
机会节点:中间过程
叶节点:最终的决策结果(决策树最后所有的输出都在叶节点)

决策树思想

决策树决策的过程,其实可以理解为你(决策树)和样本之间的对话,你不停的去询问样本的样本属性,这个询问过程,遵循的规则就是,一、你每次询问只能询问样本点的一个维度的属性,例如你可以问样本的颜色,或者大小,要想询问大小和颜色必须分两次询问,你不能同时询问样本大小和颜色.二、你每次询问的维度必须找最优的维度,不能乱问.三、在询问完成之后,样本回答了,你要按照最优的划分去让样本进入下一个对应的询问.这样一来,对于一个样本就可以不停的沿着一条询问路径走下去,落到最终的分类里面.这个过程有点像之前很火的微软小冰猜偶像的过程.至于为什么成为树,是因为,这样询问下去,会出现一个金字塔一样的树结构.

决策树算法流程

具体算法流程如下:


决策树算法流程.png

构建决策树关键

构建决策树的关键在以下两点:
1、构造最佳节点和最佳分支(也就是例子中每一次询问时,我们要找到询问最优的维度,和最优的分裂进入下一次询问)
2、决策树剪枝策略(这个是为了防止过拟合出现)。

1、构造最佳节点和最佳分支(也就是例子中每一次询问时,我们要找到询问最优的维度,和最优的分裂进入下一次询问)

对于第一个问题(其实就是算法中的第8步)。采用贪心算法,每一步都构造当前最好的节点和分支,怎么构造最好的节点和分支,常用一下指标:
a、信息增益。
b、信息增益率。
c、GINI系数。
d、误差平方和(最小二乘法)。

信息增益

信息增益的概念来自信息熵。信息熵的定义如下:
Ent(D)=- \sum_{k=1}^Kp_klog_2(p_k)

其中:D为一个集合,p_k为集合中元素的概率,k为集合中元素的类别,为什么取负号,是因为p_k<=1,所以log_2p_k<=0,从而加上负号之后,保证信息熵为正值。

由定义可以看出。信息熵是单调函数,所以,取极端值分析,若D里面的东西都是一模一样的,那么p_k=1,从而Ent(D)等于0,所以,信息熵越小,说明集合中信息的纯度越高,信息熵越大,集合中越杂乱,也可以认为信息熵的大小代表信息的不确定性,信息熵越大,信息越不确定。

有了信息熵的概念,可以像概率一样,定义条件熵,联合熵的概念。条件熵的概念来自概率的条件概率,定义如下:
条件熵 Ent(Y|X) 定义为 X 给定条件下 Y 的条件概率分布的熵对X的数学期望:
Ent(Y|X)=\sum_xp(x_i)Ent(Y|X=x_i)=-\sum_xp(x_i)\sum_yp(y_i|x_i)log_2(p(y_i|x_i))=-\sum_x\sum_yp(x_i,y_i)log_2(p(y_i|x_i))

联合熵Ent(X,Y)定义为:
Ent(X,Y)=-\sum_x\sum_yp(x_i,y_i)log_2(p(x_i,y_i))

有了信息熵的概念,就可以定义信息增益的概念,给定集合D任意一种划分a,集合D本身具有一个原始的信息熵,在知道a的前提下D就有了一个条件熵,那么给定a和没有给a的连个熵的差就是给定划分a的信息增益。

公式如下:
Gain(D,a)=Ent(D)-Ent(D|a)=Ent(D)- \sum_m \frac{|D_m|}{|D|} Ent(D_m)

其中:{D_m,m=1,2,3..M} 是划分a对集合的划分。

有了信息增益,在决策树的每一个机会节点,只需要要追求当期信息增益最大的方向去生长子节点或者叶子即可。

但是信息增益是个绝对值,会出现一个问题,就是分类数目越多时,信息增益越大,既对样本分的很细时,信息增益越大。如果样本分的过细,是没有任何意义的,是属于过拟合的,例如:一个样本集需要分成n类,如果给这个样本集随便进行一个编号(1,2,...M),如果把这个编号作为样本的一个属性,求其信息增益,那么这个信息增益会是1,非常的大。但是这个编号无意义,很明显出现了过拟合。

信息增益率

用信息增益去构建决策树,是ID3,C4.0决策树使用的,在C4.5决策树中为解决上面说的问题,没有直接使用信息增益,而是用了相对的信息增益率,其信息增益率的定义如下:
Gain\_ratio(D,a)=\frac{Gain(D,a)}{IV(a)}

其中:IV(a)=- \sum_{v=1}^V \frac{|D_v|}{|D|}log_2(\frac{|D_v|}{|D|}) 其称为a的固有值(intrinsic value),可以理解为a的信息熵,可以看出,若分类方式a分的类别越多,那么他的值也就越大,从而可以避免使用信息增益是产生偏向分类多的问题。从而避免了过拟合问题。

上面出现了三种决策树,ID3 、C4.0、C4.5,还有一种决策树叫做CART(分类回归决策树),分类回归决策树,既可以做分类问题,也可以做回归问题。用CART做分类问题,在节点判断当前分类最优时,选用基尼指数作为判断依据。用CART做回归问题,在节点判断当前分类最优时,选用误差平方和作为判断依据,

基尼指数

在定义基尼指数前,要先定义基尼值。基尼值用来衡量集和D的纯度,基尼值越大说明集和越杂乱,基尼值越小,说明集和纯度越高。其定义公式如下:
Gini(D)=\sum_{k=1}^K \sum_{k'=k}^K p_k p_{k'}=1-\sum_{k=1}^Kp_k^2

其中:K表示集合D被划分的类别.p_k是在D中第k类被选出的概率.

所以,集合的基尼值可以看作,从集合中随机抽取两个样本,其类别不一致的概率。

所以基尼指数定义如下:
Gini_index(D,a)=\sum_{v=1}^V \frac{|D^v|}{|D|}Gini(D^v)

其中:a表示一种划分集合的方式,V表示划分集合的个数,||表示集合中元素的个数。所以基尼指数可以理解为,集合D在给定一个划分a下的划分集合的基尼值得数学期望。

由基尼指数得定义可以看出,基尼指数越小,这个划分越好.如果不好理解,可以去极端值,这个划分刚好是完美划分,意思就是,每个划分得小集合D^v里面都只有一个类别,所以,Gini(D^v)=0,所以,Gini_index(D,a)=0,所以基尼指数越小,划分越好,我们要寻找基尼指数最小的划分,作为当前最优划分.

误差平方和

对于一些连续得数值,也就是非类别型变量,CART得处理方式比较特别,是使用二分法,一步一步去逼近真实值,也就是对区间不停的进行二分类,直到区间的中点和真实值的误差在运行的范围内.所以其本质还是一个分类问题.只不过最后输出不是一个类别,而是一个集体的数值.

对于每一课决策树,最后得形式都可以表示如下:
f(x)=\sum_{m=1}^MC_mI(x\in R_m)

其中:x是样本点,M是最后分成得类别个数,m是每一个类得类别编号.C_m在分类问题中表示当前类别属于R_m类别得概率,C_m在回归问题中,是最后优化出来得具体得值.I(x\in R_m)是机器学习中最长用得一种表示方式,其意思为:样本x若在R_m中则是1,若不在R_m中则是0.这种写法非常方便,可以把机器学习的问题和结论全部公式化,推导起来非常方便,但一个最大的问题是不易理解.

有了上面决策树最后的函数,那么对于回归问题,就可以定义误差,其误差表达式为:
J=\sum_i^N(y_i-f(x_i))^2=\sum_i^N(y_i-\sum_{m=1}^MC_mI(x_i\in R_m))^2

我们从上面的公式反推出C_m的具体值.我们要获得最优的回归树,优化方法就是要是上面的公式最小,即最小化误差平方和(最小二乘法):
min(J)=min\sum_i^N(y_i-f(x_i))^2=min\sum_i^N(y_i-\sum_{m=1}^MC_mI(x_i\in R_m))^2

因为C_m仅仅和R_m相关,即一个类别R_m对应一个值C_m,所以,要使上式最小,对C_m求导得到:
\frac{\partial J}{\partial C_m}=\sum_{i=1}^N 2(y_i-\sum_{m=1}^MC_mI(x_i\in R_m))I(x_i\in R_m)

对于固定的m_0,有
\frac{\partial J}{\partial C_{m_0}}=\sum_{k=1}^K 2(y_i-\sum_{k=1}^KC_{m_0})

所以,J取极小值时,C_{m_0}=\frac{\sum_{k=i}^Ky_i}{K}

其中:K为类别R_{m_0}中的样本数量,y_i为样本点x_i对应的真实值.

所以,我们在决策树的每一个决策节点,只需要找到的一种区间划分,满足划分后的每一个区间的中点和真实值的差的平方和最小既可以.

即每次树节点分裂根据以下公式的j,s
\min_{j,s}[\min \sum_{x_i \in R_1(j,s)}(y_i-C_1)^2 +\min \sum_{x_i \in R_2(j,s)}(y_i-C_2)^2]

其中:j是样本点的第j个属性,s是样本点在某个属性上的切分点,R_1(j,s)是样本被切开的第一个集合,切的规则是按第j个属性,用s为切点切割.R_2(j,s)是样本点被刚刚切剩下的集合.C_1=avg(y_i) 此时x_i\in R_1(j,s) , C_2=avg(y_i) 此时x_i\in R_2(j,s)

至此,我们解决了构建决策树时第一个关键.解决了这个问题,我们就可以构建一棵决策树.但是棵决策树很容易过拟合,为解决过拟合问题,我们引入下一个问题,决策树剪枝策略.

2、决策树剪枝策略(这个是为了防止过拟合出现)。

决策树出现过拟合,是因为分的过细,最后导致产生的叶子太多了,每个叶子上的样本比例过少.所以我们应该控制决策树产生过多的叶节点.也就是对决策树进行剪枝.

决策树的剪枝策略分为两种,一种时预剪枝,一种是后剪枝.顾名思义,预剪枝是在训练决策树时,在机会节点是否向下分裂时,添加一步判断,若满足这个条件就不进行分裂,直接成为一个叶节点.若不满足正常进行.而后剪枝是,在决策树训练完成后,对已形成的决策树自下而上的每个节点进行逐一裁剪,添加判断条件,判断剪掉该节点是变好了还是变坏了,若变好,则剪枝.

所以不难看,预剪枝减少了计算量和计算时间,但是由于每一步追求的是当前最优,所以预剪枝会出现欠拟合的问题.后剪枝计算时间和计算量都比较大,但效果会比预剪枝好.不论是预剪枝还是后剪枝,其关键在于判断条件.

因为是解决过拟合问题,这个判断条件可直接选择在验证集上的精度.若验证集剪枝后精度提高,则剪枝.精度下降则步剪枝.精度的定义如下:
acc(f;D)=\frac{1}{m}\sum_{i=1}^mI(f(x_i=y_i))

其中:f为决策树,D为样本集合.m是岩本集合D的样本个数.

实例

下面是对一个基金数据集的决策树的例子.

from itertools import product
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from sklearn import datasets
from sklearn.tree import DecisionTreeClassifier
 

# 仍然使用自带的iris数据
#iris = datasets.load_iris()

iris = pd.read_csv(r'fund.txt',header=None)
iris.columns =  ['Classlabel', 'FundScale', 'FundSharpeRatio']


#X = iris.data[:, [0, 1]]

#y = iris.target
X=iris[['FundScale', 'FundSharpeRatio']]
y=iris['Classlabel']
#print(iris.feature_names)
# 训练模型,限制树的最大深度4
clf = DecisionTreeClassifier(max_depth=4)
#拟合模型
clf.fit(X, y)


# 画图
x_min, x_max =np.array(X['FundScale']).min() - 1, np.array(X['FundScale']).max() + 1
y_min, y_max = np.array(X['FundSharpeRatio']).min() - 1, np.array(X[ 'FundSharpeRatio']).max() + 1

print(x_min, x_max)
print(y_min, y_max)
xx, yy = np.meshgrid(np.arange(x_min, x_max, 200000),np.arange(y_min, y_max, 1))

Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

plt.contourf(xx, yy, Z, alpha=0.4)
plt.scatter(np.array(X['FundScale']), np.array(X[ 'FundSharpeRatio']), c=y, alpha=0.8)
plt.show()
29274.988 15265537319.592
-1.5773000000000001 29.8519
output_59_1.png
DecisionTreeClassifier()函数为创建一个决策树模型,其函数的参数含义如下所示:

·criterion:gini或者entropy,前者是基尼系数,后者是信息熵。 所以这个参数就可以指定是ID3或者C4.5决策树还是CRAT决策树,但是没有C4.5决策树.
·splitter: best or random 前者是在所有特征中找最好的切分点 后者是在部分特征中,默认的”best”适合样本量不大的时候,而如果样本数据量非常大,此时决策树构建推荐”random” 。
·max_features:None(所有),log2,sqrt,N 特征小于50的时候一般使用所有的
·max_depth: int or None, optional (default=None) 设置决策随机森林中的决策树的最大深度,深度越大,越容易过拟合,推荐树的深度为:5-20之间。
·min_samples_split:设置结点的最小样本数量,当样本数量可能小于此值时,结点将不会在划分。
·min_samples_leaf: 这个值限制了叶子节点最少的样本数,如果某叶子节点数目小于样本数,则会和兄弟节点一起被剪枝。
·min_weight_fraction_leaf: 这个值限制了叶子节点所有样本权重和的最小值,如果小于这个值,则会和兄弟节点一起被剪枝默认是0,就是不考虑权重问题。
·max_leaf_nodes: 通过限制最大叶子节点数,可以防止过拟合,默认是"None”,即不限制最大的叶子节点数。
·class_weight: 指定样本各类别的的权重,主要是为了防止训练集某些类别的样本过多导致训练的决策树过于偏向这些类别。这里可以自己指定各个样本的权重,如果使用“balanced”,则算法会自己计算权重,样本量少的类别所对应的样本权重会高。
·min_impurity_split: 这个值限制了决策树的增长,如果某节点的不纯度(基尼系数,信息增益,均方差,绝对差)小于这个阈值则该节点不再生成子节点。即为叶子节点 。

#正确安装graphviz
from IPython.display import Image  
from sklearn import tree
import pydotplus 
dot_data = tree.export_graphviz(clf,out_file=None,feature_names=[ 'FundScale', 'FundSharpeRatio'],class_names='Classlabel',filled=True, rounded=True,special_characters=True)  
graph = pydotplus.graph_from_dot_data(dot_data)  
Image(graph.create_png())
output_61_0.png

至此决策树的介绍就结束了.下面介绍booststrap.

booststrap(自助法)

booststrap在统计学中是一种统计方法,其具体操作为,对于一个样本为N的集合,对这个集合进行有放回的抽样实验,抽出M个样本,这M个样本形成新的统计样本集,对于新的样本集会有一些列统计指标,用这些统计指标代替原始原本的统计指标.

booststrap在机器学习中狭义的指,对于机器学习的训练集和测试集的一种验证划分方法.在机器学习中训练集和测试集有三种划分方法:留出法,k折交叉验证法,还有booststrap(自助法).每种方法都有优缺点.

留出法

这种方法简单直接,是很生硬的将数据集人为的分为训练集和测试集两个不相交的集合.这种方法优点是简单直接,但的问题在于训练样本和测试样本数量权衡问题,训练样本过多,必然导致测试样本减少,从而导致结果的方差较大,结果不够稳定,如果测试样本过多,可以保证结果的稳定性,但是训练集样本就会过少,导致结果的偏差较大.所以如何衡量训练集和测试集样本数量是这种方法的主要问题.常用的做法是训练集保证全部样本集的2/3~4/5.

k重交叉验证

这种方式实际上实对上面方式的一种弥补,首先对样本集分成k个互斥的集合,这里要保证每个集合的分布尽可能一致,然后在这k 个互斥的集合中选择一个作为测试集,其他的合在一起作为训练集,以此训练一个模型,然后再重复这个步骤,再k个集合中选择一个刚才没有选择的集合作为测试集,其他的合在一起作为训练集.再重新训练一个模型出来,这样既可以训练k个模型出来,取这k个模型的均值.作为最后的结果.这种方式解决了偏差问题.但是缺点为计算量会非常大.

k重交叉验证中有一个特殊情况就是留一法,意思是保留一个样本作为训练集,其他的作为测试集,多次训练,让每个样本都有做测试集的机会.可想而知这种方法虽然准确稳定,但是计算量超级大.

自助法(booststrap)

自助法是在原数据集中进行随机采样,随机采样得到的样本集为训练集,剩下的为测试集.自助法有个问题就是随机选择的样本和原始样本的分布可能不一样,所以也加大了偏差.且自助法只适用于样本数量少的训练.

对于booststrap,每一个样本不被选中的概率为以下:
P=(1-\frac{1}{m})^m
其中:m是随机选出样本的个数。

这个值的极限(m趋于无穷大时)为\frac{1}{e}=36.8\%,所以有大约36.8%的样本没有被随机选出来。

Bagging(袋装法)

这种方法和boosting(提升法)相对应,两者都是一种集成学习的方法.集成学习,就是把弱分类器放在一起搞出一个强分类器的学习过程.

bagging的思想是,重复多次对原始样本进行boostrap(自助法)训练模型,得到多个模型,若是分类问题,用这些模型进行投票,若是回归问题,用这些模型结果求均值.

bagging是多次对原样本集进行自助抽样训练模型,每次自助抽样产生的样本集为训练集,剩下的为测试集,比例大约是63:37。所以每次的训练集既有公共本分又要不一样的部分,所以训练出的模型也不一样,所以可以保证泛化性,防止过拟合的出现。

如上面所说,每次训练约有1/3的样本没有参与模型训练,但这1/3的样本还有其他用处,这1/3的样本,被称为袋外样本(Out-Of-Bagging),这些样本并非无用,这些样本一般来求袋外数据误差(也称为Bagging泛化误差),这个误差有什么用呢,这个误差可直接用来判断模型的泛化性,不需要额外的测试集进行测试。

袋外数据误差(OOB error)定义如下:
Error^{oob}=\frac{1}{D}\sum_{(x,y)\in D}I(H^{oob}(x)\neq y)

其中D为原始样本集,
H^{oob}(x)=argmax_{y \in Y }\sum_{t=1}^TI(h_t(x)=y)I(x \notin D_t) , 实际为在袋外样本的bagging方法的最后预测结果,
t为第t次使用自助抽样训练模型,
y为label,
T 训练模型的次数,
h_t为第t个预测模型,
D_t为第t次自助抽样时,抽样出的样本集。

袋外泛化误差(OOB error)的本质其实就是袋外样本在对应训练模型上分错的数量除以原始样本的总数量。

注:OOB Error是随机森林泛化误差的一个无偏估计,它的结果近似于需要大量计算的k折交叉验证。

bagging的算法流程如下
bagging.jpg

其中:D_{bs}为自助随机抽样产生的样本集。

boosting(提升法)

boosting的思想是,每次训练都是用全部样本,只不过下次训练用的样本的权重会被改变,具体的做法为,第一次每个样本的权重是等权,以此训练出第一个模型,得到第一个模型后,根据每个样本模型预测的值和真实值之间的误差来调节各个样本的权重,误差越大的样本权重越大,以此循环,直到训练完成指定的步数。最后,将各个步骤训练的模型做一个加权和。boosting的代表算法是Adaboost,所以具体的算法步骤和推导在Adaboost中介绍。

boosting和bagging都是集成学习的方法,但两者的方向不一样.boosting主要是解决偏差问题,偏差很小,但是方差很大,bagging主要是解决方差问题,方差很小,但是偏差很大.而且,我们可以看出,bagging是一种并行的方式,模型和模型间没有依赖性,彼此独立.boosting是一种串行的方式,模型和模型间是一种强的依赖关系.bagging的代表是随机森林,boosting的代表是Adaboost和GBDT.

既然讲到偏差和方差,那这里再给出偏差和方差的概念,以及偏差方差分解理论,这个概念在机器学习中非常重要。

偏差方差分解理论
方差

var(x)=E_D[(f(x:D)-\overline f(x))^2]

噪音

e^2=E[(y_D-y)^2]

偏差

bias^2(x)=(\overline f(x)-y)^2

其中:D为训练样本点集合,x为样本,y_D为样本的标签,y为样本的真实标签。y_Dy的区别是样本测量记录值和真实情况值的区别,例如一个物体长1.23456789m,但是测量的尺子精确度到cm,所以在记录样本时,就会写1.235m.从而产生了误差。这个误差的期望是零,而且和样本、模型完全没有关系。
\overline f(x)=E_D[f(x;D)],表示模型在多个训练集(样本数一样的不同训练集)上的预测值的均值,注意这个地方不是对所有的x取均值,而是对不同的训练集训练出得模型在同一个点上取均值。

由上面方差、噪音和偏差的定义,我们可以看出方差来衡量模型的稳定性,偏差是用来衡量模型的准确性,噪音来衡量模型学习的难度。两个指标再加上噪音,就可以定义模型的泛化误差。

偏差和方差的直观理解,可以看下图

[图片上传失败...(image-909cc4-1578319182808)]

泛化误差

E(f;D)=E_D[(f(x;D)-y_D)^2]=E_D[(f(x;D)- \overline f(x)+ \overline f(x)-y_D)^2]=E_D[((f(x;D)- \overline f(x))^2+(\overline f(x)-y_D)^2]-2E_D[(f(x;D)- \overline f(x))(\overline f(x)-y_D)]

其中对于E_D[(f(x;D)- \overline f(x))(\overline f(x)-y_D)]

因为\overline f(x)是一个常数,且y_D是实际测量标签,所以与f(x;D)独立,所以(\overline f(x)-y_D)(f(x;D)- \overline f(x))
独立。由独立的性质:
E_D[(f(x;D)- \overline f(x))(\overline f(x)-y_D)] =E_D[(f(x;D)- \overline f(x))]E[((\overline f(x)-y_D))]=[E_D(f(x;D)- E_D\overline f(x))]E[((\overline f(x)-y_D))]=[E_D(f(x;D)- \overline f(x))]E[((\overline f(x)-y_D))]=0

所以:
E(f;D)=E_D[(f(x;D)- \overline f(x))^2] + E_D[(\overline f(x)-y_D)^2]

所以继续推导:
E(f;D)=E_D[(f(x;D)- \overline f(x))^2] + E_D[(\overline f(x)- y+y- y_D)^2] =E_D[(f(x;D)- \overline f(x))^2]+E_D[(\overline f(x)- y)^2]+E_D[(y- y_D)^2]-2E_D [(\overline f(x)- y)(y- y_D)]

其中对于2E_D[(\overline f(x)- y)(y- y_D)]
因为\overline f(x)是常数,y-y_D是测量误差活着噪音,与y本身没有任何关系,所以y-y_D\overline f(x)- y相互独立,所以有
E_D[(\overline f(x)- y)(y- y_D)]=E_D[(\overline f(x)- y)]E_D[(y- y_D)],而噪音的期望为零,所以这一项总体也为零。

所以泛化误差为:
E(f;D) =E_D[(f(x;D)- \overline f(x))^2]+E_D[(\overline f(x)- y)^2]+E_D[(y- y_D)^2] =var(x)+bias^2(x)+e^2

所以得出结论泛化误差可以分解为偏差,方差和噪音的和。泛化误差可以综合考虑模型的稳定性和准确性,所以一般直接用泛化误差就可以判断模型的好坏。

对于训练一个模型,我们用训练强度来指,模型对训练集的拟合程度,如过,模型的训练强度太低,模型欠拟合,这样会导致模型的偏差过大,而方差过小,此时,偏差主导了泛化误差,随着训练程度的加大,模型拟合越来越好,偏差会越来越小,但是方差会越来越大,当训练强度过大时,模型会出现过拟合,此时,方差主导了泛化误差。

具体泛化误差,偏差,方差关系如下图所示:

[图片上传失败...(image-74d319-1578319182808)]

上面是插入的一段关于机器学习方差偏差分解的理论,下面继续接着Bagging介绍,有了bagging,自然而然的就会引出随机森林的算法.下面介绍随机森林的算法.

随机森林

随机森林,顾名思义两个重点,随机和森林,随机主要是两个地方随机,森林是说里面有一堆决策树.应用Bagging的方法,对原始样本集,随机选取M个样本,再再这M个样本上随机选取N个属性(连个随机),形成一个新的数据集,用这个新的数据集训练一颗决策树,并且重复这个步骤K次,那么一个简单的随机深林生成了.这个随机森林中有K棵决策树.

由上面两个随机的抽样可以看出.随机森林缓解了树模型方差大的特点,方差减小,而且方差的减小补偿了偏差的增大,因此总体而言是更好的模型。
在构建决策树的时候,.随机森林的每棵决策树都最大可能的进行生长而不进行剪枝;在对预测输出进行结合时,.随机森林通常对分类问题使用简单投票法,回归任务使用简单平均法。 由于随机森林使用的是baaging算法,所以不需要交叉验证,OOB Error(袋外数据误差)是随机森林泛化误差的一个无偏估计,它的结果近似于需要大量计算的k折交叉验证。

优点

1)在数据集上表现良好,相对于其他算法有较大的优势(训练速度、预测准确度)。
2)能够处理很高维的数据,并且不用特征选择,而且在训练完后,给出特征的重要性。
3)容易做成并行化方法。

缺点

1)在噪声较大的分类或者回归问题上会过拟合。
2)离群点敏感性。

噪声点影响:Adaboost > SVM > RF,IForest(孤立森林)可以做异常检测。

随机森林的流程图
随机森林.png

流程说明:
1、从原始训练集中使用Bootstraping方法随机有放回采样取出m个样本,共进行n_tree次采样。生成n_tree个训练集.
2、对n_tree个训练集,我们分别训练n_tree个决策树模型.
3、对于单个决策树模型,假设训练样本特征的个数为n,那么每次分裂时根据信息增益/信息增益比/基尼指数 选择最好的特征进行分裂.
4、每棵树都这样分裂下去,直到该节点的所有训练样例都属于同一类。在决策树的分裂过程中不需要剪枝.
5、将生成的多颗决策树组成随机森林。对于分类问题,按照多棵树分类器投票决定最终分类结果;对于回归问题,由多颗树预测值的均值决定最终预测结果.
注意:OOB(out-of-bag ):每棵决策树的生成都需要自助采样,这时就有1/3的数据未被选中,这部分数据就称为袋外数据。用这些数据根据袋外数据误差公式计算袋外数据误差,用袋外数据误差评估随机森林效果.

特征选择

随机森林除了可以进行分类和回归问题意外,还有一个非常重要的作用,就是特征选择.随机森林,可以在训练的过程中,根据特征的重要性,对特征进行排序,选出重要的特征.

用随机森林进行特征重要性评估的思想就是看每个特征在随机森林中的每棵树上做了多大的贡献,然后取个平均值,最后比一比特征之间的贡献大小。
贡献大小通常使用基尼指数(Gini index)或者袋外数据(OOB)错误率作为评估指标来衡量。
这里我们介绍一下基尼指数来评价的方法。
我们将变量重要性评分(variable importance measures)用VIM来表示,将Gini指数用GI来表示,假设m个特征a_1,a_2,a_3,......a_m,现在要计算出每个特征a_j的Gini指数评分VIM_j(Gini),即第j个特征在随机森林中所有决策树中节点分裂不纯度的平均改变量。

第一步:我们找到随机森林中所有决策树上以a_j为特征进行分裂的节点,把这些节点放在集合M中,对于M中的每一个节点,分别计算它与由它分裂的叶节点的基尼值之差,例如:对于节点i,计算此节点的基尼值
Gini(i)=1-\sum_{k=1}^Kp_k
其中:K是节点在特征a_j上能划分出的种类,也就是此节点分裂出叶节点个数.然后计算每一个叶节点的Gini指数.记作Gini(i)_k ,(k=1,2,3...K)
所以节点i最终的基尼值之差为:
Gini(i)-\sum_{k=1}^KGini(i)_k

第二步:把集合M中的每一个节点的基尼值之差全部按第一步的公式计算出来.求和得到特征a_j的重要性.即:
VIM_j=\sum_{i=1}^I( Gini(i)-\sum_{k=1}^KGini(i)_k)

第三步:用上述两步的方法求得每一个特征的重要性,让后进行归一化,最后降序排序,

下面是一个简单的用随机森林确定特征重要性的列子.

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier

#url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data'
url1 = pd.read_csv(r'fund2.txt',header=None)
# url1 = pd.DataFrame(url1)
# df = pd.read_csv(url1,header=None)
url1.columns =  ['Class label', 'FundScale', 'FundSharpeRatio', 'FundManagerMonth',
              'FundMonth', 'FundDownsideDev', 'FundOmegaRatio',
              'FundOmegaRatioIntegral', 'FundUpCapture', 'FundDownCapture']

# 查看几个标签
# Class_label = np.unique(url1['Class label'])
# print(Class_label)
# 查看数据信息
# info_url = url1.info()
# print(info_url)

# 除去标签之外,共有13个特征,数据集的大小为178,
# 下面将数据集分为训练集和测试集
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import Imputer

# url1 = url1.values
# x = url1[:,0]
# y = url1[:,1:]
x,y = url1.iloc[:,1:].values,url1.iloc[:,0].values
x_train,x_test,y_train,y_test = train_test_split(x,y,test_size=0.3,random_state=0)
feat_labels = url1.columns[1:]

# n_estimators:森林中树的数量
# n_jobs  整数 可选(默认=1) 适合和预测并行运行的作业数,如果为-1,则将作业数设置为核心数
forest = RandomForestClassifier(n_estimators=10000, random_state=0, n_jobs=-1)

x_train = Imputer().fit_transform(x_train)

 
forest.fit(x_train, y_train)

# 下面对训练好的随机森林,完成重要性评估
# feature_importances_  可以调取关于特征重要程度
importances = forest.feature_importances_
print("重要性:",importances)
x_columns = url1.columns[1:]
indices = np.argsort(importances)[::-1]
for f in range(x_train.shape[1]):
    # 对于最后需要逆序排序,我认为是做了类似决策树回溯的取值,从叶子收敛
    # 到根,根部重要程度高于叶子。
    print("%2d) %-*s %f" % (f + 1, 30, feat_labels[indices[f]], importances[indices[f]]))


# 筛选变量(选择重要性比较高的变量)
threshold = 0.15
x_selected = x_train[:,importances > threshold]
# 可视化
import matplotlib.pyplot as plt
plt.figure(figsize=(10,6))
plt.title("Importance of features",fontsize = 18)
plt.ylabel("import level",fontsize = 15,rotation=90)
plt.rcParams['font.sans-serif'] = ["DejaVu Sans"]
plt.rcParams['axes.unicode_minus'] = False
for i in range(x_columns.shape[0]):
    plt.bar(i,importances[indices[i]],color='orange',align='center')
    plt.xticks(np.arange(x_columns.shape[0]),x_columns[indices],rotation=90,fontsize=15)
plt.show()
重要性: [0.14694769 0.10870715 0.11751638 0.07467818 0.10645336 0.11860345
 0.09992749 0.10804669 0.1191196 ]
 1) FundScale                      0.146948
 2) FundDownCapture                0.119120
 3) FundOmegaRatio                 0.118603
 4) FundManagerMonth               0.117516
 5) FundSharpeRatio                0.108707
 6) FundUpCapture                  0.108047
 7) FundDownsideDev                0.106453
 8) FundOmegaRatioIntegral         0.099927
 9) FundMonth                      0.074678
output_128_1.png

下面给出一个随机森林在分类算法上的例子

#_*_coding:UTF_8_*_

# 导入需要导入的库
import pandas as pd
import  numpy as np
import math
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn import model_selection ,metrics
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.datasets import make_blobs
import warnings

# 忽略一些版本不兼容等警告
warnings.filterwarnings("ignore")

# 每个样本有几个属性或者特征
n_features = 2
#x,y = make_blobs(n_samples=300,n_features=n_features,centers=6)
x,y = url1.iloc[:,1:3].values,url1.iloc[:,0].values
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state=1,train_size=0.7)

# 绘制样本显示
# plt.scatter(x[:,0],x[:,1],c=y)
# plt.show()

# 传统决策树,随机森林算法 极端随机数的区别
DT = DecisionTreeClassifier(max_depth=None,min_samples_split=2,random_state=0)

RF = RandomForestClassifier(n_estimators=10,max_features=math.sqrt(n_features),
                            max_depth=None,min_samples_split=2,bootstrap=True)

EC = ExtraTreesClassifier(n_estimators=10,max_features=math.sqrt(n_features),
                          max_depth=None,min_samples_split=2,bootstrap=False)



# 训练
DT.fit(x_train,y_train)
RF.fit(x_train,y_train)
EC.fit(x_train,y_train)

#区域预测
# 第0列的范围
x1_min,x1_max = x[:,0].min(),x[:,0].max()
# 第1列的范围
x2_min,x2_max = x[:,1].min(),x[:,1].max()
# 生成网格采样点行列均为200点
x1,x2 = np.mgrid[x1_min:x1_max:200j,x2_min:x2_max:200j]
# 将区域划分为一系列测试点用去学习的模型预测,进而根据预测结果画区域
area_sample_point = np.stack((x1.flat,x2.flat),axis=1)


# 所有区域点进行预测
area1_predict = DT.predict(area_sample_point)
area1_predict = area1_predict.reshape(x1.shape)


area2_predict = RF.predict(area_sample_point)
area2_predict = area2_predict.reshape(x1.shape)

area3_predict = EC.predict(area_sample_point)
area3_predict = area3_predict.reshape(x1.shape)


# 用来正常显示中文标签
mpl.rcParams['font.sans-serif'] = [u'DejaVu Sans']
# 用来正常显示负号
mpl.rcParams['axes.unicode_minus'] = False

# 区域颜色
classifier_area_color = mpl.colors.ListedColormap(['#A0FFA0','#FFA0A0','#A0A0FF'])
# 样本所属类别颜色
cm_dark = mpl.colors.ListedColormap(['r','g','b'])

# 绘图
# 第一个子图
plt.subplot(2,2,1)

plt.pcolormesh(x1,x2,area1_predict,cmap = classifier_area_color)
plt.scatter(x_train[:,0],x_train[:,1],c =y_train,marker='o',s=50,cmap=cm_dark)
plt.scatter(x_test[:,0],x_test[:,1],c =y_test,marker='x',s=50,cmap=cm_dark)

plt.xlabel('data_x',fontsize=8)
plt.ylabel('data_y',fontsize=8)
plt.xlim(x1_min,x1_max)
plt.ylim(x2_min,x2_max)
plt.title(u'DecisionTreeClassifier',fontsize=8)
plt.text(x1_max-9,x2_max-2,u'o-------train ; x--------test$')

# 第二个子图
plt.subplot(2,2,2)

plt.pcolormesh(x1,x2,area2_predict,cmap = classifier_area_color)
plt.scatter(x_train[:,0],x_train[:,1],c =y_train,marker='o',s=50,cmap=cm_dark)
plt.scatter(x_test[:,0],x_test[:,1],c =y_test,marker='x',s=50,cmap=cm_dark)

plt.xlabel('data_x',fontsize=8)
plt.ylabel('data_y',fontsize=8)
plt.xlim(x1_min,x1_max)
plt.ylim(x2_min,x2_max)
plt.title(u'RandomForestClassifier',fontsize=8)
plt.text(x1_max-9,x2_max-2,u'o-------train ; x--------test$')

# 第三个子图
plt.subplot(2,2,3)

plt.pcolormesh(x1,x2,area3_predict,cmap = classifier_area_color)
plt.scatter(x_train[:,0],x_train[:,1],c =y_train,marker='o',s=50,cmap=cm_dark)
plt.scatter(x_test[:,0],x_test[:,1],c =y_test,marker='x',s=50,cmap=cm_dark)

plt.xlabel('data_x',fontsize=8)
plt.ylabel('data_y',fontsize=8)
plt.xlim(x1_min,x1_max)
plt.ylim(x2_min,x2_max)
plt.title(u'ExtraTreesClassifier',fontsize=8)
plt.text(x1_max-9,x2_max-2,u'o-------train ; x--------test$')

# 第四个子图
plt.subplot(2,2,4)
y = []
# 交叉验证
score_DT = cross_val_score(DT,x_train,y_train)
y.append(score_DT.mean())
score_RF = cross_val_score(RF,x_train,y_train)
y.append(score_RF.mean())
score_EC = cross_val_score(EC,x_train,y_train)
y.append(score_EC.mean())

print('DecisionTreeClassifier交叉验证准确率为:'+str(score_DT.mean()))
print('RandomForestClassifier交叉验证准确率为:'+str(score_RF.mean()))
print('ExtraTreesClassifier交叉验证准确率为:'+str(score_EC.mean()))

x = [0,1,2]
plt.bar(x,y,0.4,color='green')
plt.xlabel("0--DecisionTreeClassifier;1--RandomForestClassifier;2--ExtraTreesClassifie", fontsize=8)
plt.ylabel("avg_acc", fontsize=8)
plt.ylim(0.7, 0.99)
plt.title("cross", fontsize=8)
for a, b in zip(x, y):
    plt.text(a, b, b, ha='center', va='bottom', fontsize=10)
plt.show()
DecisionTreeClassifier交叉验证准确率为:0.7582268679829656
RandomForestClassifier交叉验证准确率为:0.7663569492837786
ExtraTreesClassifier交叉验证准确率为:0.7423538521099496
output_130_1.png

所以从上面的两个例子中,可以看出,实现一个随机森林,最重要的是,在随机森林训练模型时对其调参.

随机森林调参

n_estimators=10:决策树的个数,越多越好,但是性能就会越差,至少100左右(具体数字忘记从哪里来的了)可以达到可接受的性能和误差率。

bootstrap=True:是否有放回的采样。

oob_score=False:oob(out of band,带外)数据,即:在某次决策树训练中没有被bootstrap选中的数据。多单个模型的参数训练,我们知道可以用cross validation(cv)来进行,但是特别消耗时间,而且对于随机森林这种情况也没有大的必要,所以就用这个数据对决策树模型进行验证,算是一个简单的交叉验证。性能消耗小,但是效果不错。

n_jobs=1:并行job个数。这个在ensemble算法中非常重要,尤其是bagging(而非boosting,因为boosting的每次迭代之间有影响,所以很难进行并行化),因为可以并行从而提高性能。1=不并行;n:n个并行;-1:CPU有多少core,就启动多少job

warm_start=False:热启动,决定是否使用上次调用该类的结果然后增加新的。

class_weight=None:各个label的权重。

Adaboost

Booting的第一个代表算法是Adaboost,其主要思想是,全样本多次串行训练,每次训练时,虽然训练样本集一样,但是其每个样本的权重不一样,每次训练样本的权重根据上步训练结果得到,上步训练结果的损失函数大,那么权重就大,下次训练时就会重点考虑这些样本的影响。这样,串行达到指定的步数之后,将每一步训练得到的结果加权相加,就得到了最后的模型。

其流程图如下

boosting.PNG

如上图最后模型为:
F(x)= \sum_{t=1}^T(\alpha_t f_t(x))

所以可以看出模型的关键为:1、每一步训练时确定样本的权重,2、最后加总时确定每一步f_i(基模型)的权重。

对于boosting方法,因为模型是每一步串行运算,分类问题和回归问题在每一步之间推导完全不一样,所以分类模型和回归模型需要分开推导,而随机森林则不需要,下面分别给出模型的详细推导,并解决上面两个关键问题。

基于分类的Adaboost算法

机器学习中每一个模型的推导,其出发点,基本都是最小化损失函数(所以机器学习的问题其本质都是最优化)。

对于分类问题,这里考虑二分类问题,对于原始的Adaboost只支持二分类,对于多分类问题,Adaboost可以使用One-Vs-One或者One-Vs-Rest的方法转化成二分类问题来解决。

考虑二分类问题,训练集D=\{(x_1,y_1),(x_2,y_2)...(x_m,y_m)\},其中y_i\in \{-1,+1 \}D_t为Adaboost第t步训练分类器时,以t-1步的结果导致权重变化后的样本分布, 设f_t(x)为样本x在样本分布D_t上的弱分类器的预测结果。这里一定要注意,这个f_t(x)并不是t步之后的样本x的最终预测值,而是第t个分类器在样本分布D_t上的预测值,而样本x的最终预测值为:
H(x)=sign(F_t(x))=sign(\sum_{i=1}^{t}\alpha_i f_i(x))

在这个地方,周志华老师的机器学习一书,混淆了H(x)与F(x),网上很多推导认为Adaboost的最终分类器函数H(x)是加法模型,其实是错误的,H(x)并不是加法模型,而上式中的F(x)才是加法模型。这一点区分开来可以为下面的权重推导做铺垫,周志华老师在未区分H(x)与F(x)的区别,导致推导时混淆不清,所以一定区分开,F(x)和H(x)的区别。H(x)是Adaboost的最终输出结果,在推导过程和算法过程中用不到。

推导Adaboost时,损失函数使用指数损失函数,至于为什么可以使用指数损失函数,将在损失函数的讨论中专门讲解,至于为什么使用指数损失函数,而不是用其他损失函数,是因为Adaboost的关键是确定每一步训练时样本的权重和最后加总时,每一步结果的权重。只有使用指数损失函数进行推导,这两个权重才能得以确定。

指数损失函数如下:
L(y,f(x))=e^{-yf(x)}
其中y为真实值,f(x)为分类器预测值。

因为Adaboost是加法模型,采用贪心算法,只追求当前一步最优,所以在当前一步时,只追求当前的期望损失最小化。即:
l_{exp}(\alpha_{t}f_t|D_t)=E_{x~D_{t}} (e^{-y \alpha_{t}f_t(x) })
对于f_t(x),是第t步对应的弱分类器的预测值,所以只能取+1,—1。这点非常重要。

对于这个期望,把样本D_t中的点分为预测对的,和预测错的,则上式可展开,为:
l_{exp}(\alpha_{t}f_t|D_t)=E_{x~D_{t}}  (e^{-y \alpha_{t}f_t(x) })=p(y=f_t(x))e^{-y \alpha_{t}f_t(x) }+ p(y \neq f_t(x))e^{-y \alpha_{t}f_t(x) } = p(y=f_t(x))e^{-\alpha_{t} }+ (1-p(y=f_t(x)))e^{ \alpha_{t} }

要是这个最小,求导,令导数等于零,得到
\alpha_t =\frac{1}{2}ln\frac{p(y=f_t(x))}{1-p(y=f_t(x))}

所以这样就解决了上面两个问题中的一个,再看在训练t步之前,样本权重的变化,即样本分布的变化。

在弱分类器权重的推导过程中,习惯用错误率来表示,上式是正确率表示的,所以上式还可以表示为:

\alpha_t =\frac{1}{2}ln\frac{1-p(y \neq f_t(x))}{p(y\neq f_t(x))}

由上式可以看出p(y \neq f_t(x))<0.5才有\alpha_t>0,所以当计算出的p(y \neq f_t(x))>0.5时,分类器可以直接放弃。

对于第t步,追求第t步的期望损失函数最小得到了t步分类器的权重,我们同时追求前t步的损失函数最小,就可推导出每步训练之前的样本分布。

对于前t步的期望损失函数(这里推导一定注意不能使用H(x)):
l_{exp}(F_t|D)=E_{x~D}  (e^{-yF_t(x) })=E_{x~D} (e^{-y (F_{t-1}(x)+\alpha_tf_t(x) ) }) =E_{x~D}  (e^{-y F_{t-1}(x)}e^{(\alpha_tyf_t(x) ) }) =e^{\alpha_t}E_{x~D}  (e^{-y F_{t-1}(x)}e^{(yf_t(x) ) })

对于e^x在零点进行二阶泰勒展开,有e^x\approx1+x+\frac{x^2}{2}

所以
l_{exp}(F_t|D) =e^{\alpha_t}E_{x~D}  (e^{-y F_{t-1}(x)}e^{(yf_t(x) ) }) =e^{\alpha_t}E_{x~D}  (e^{-y F_{t-1}(x)}(1-yf_t(x)+\frac{y^2f_t^2(x)}{2} ))

y^2=1,f_t^2(x)=1,所以:
l_{exp}(F_t|D) =e^{\alpha_t}E_{x~D}  (e^{-y F_{t-1}(x)}(1-yf_t(x)+\frac{y^2f_t^2(x)}{2} )) =e^{\alpha_t}E_{x~D}  (e^{-y F_{t-1}(x)}(1-yf_t(x)+\frac{1}{2} ))

这里我们只研究f_t(x),寻找f_t(x)使l_{exp}(F_t|D)最小。可以看出使得l_{exp}(F_t|D)最小的f_t(x) ,和使得E_{x~D} (e^{-y F_{t-1}(x)}(yf_t(x) ))最大的f_t(x)是同一个。所以,只需研究E_{x~D} (e^{-y F_{t-1}(x)}(yf_t(x) ))

即:
\underset {f_t}{argmax} E_{x~D}  (e^{-y F_{t-1}(x)}(yf_t(x) ))

又因为E_{x~D} ( e^{-y F_{t-1}(x)})是个常数。
所以:\underset {f_t}{argmax} E_{x~D}  (e^{-y F_{t-1} (x)} (yf_t(x) )) = \underset {f_t}{argmax} E_{x~D}  (\frac{ e^{-y F_{t-1} (x)}}{ E_{x~D} ( e^{-y F_{t-1}(x)}) }(yf_t(x) ))

此事对于每一个样本x,这个公式\frac{ e^{-y F_{t-1} (x)}}{ E_{x~D}  ( e^{-y F_{t-1}  (x)}  ) }是一个数值,用这个数值乘以样本的初始权重(初是分布),就会得到一个新的样本的分布,这是因为
E_{x~D}  \frac{ e^{-y F_{t-1} (x)}}{ E_{x~D} ( e^{-y F_{t-1}(x)}) }=1

把这个新的样本分布记为D_t,
既:
D_t(x)= \frac{D(x) e^{-y F_{t-1} (x)}}{ E_{x~D} ( e^{-y F_{t-1}(x)}) }

\begin{align} &\underset {f_t}{argmax} E_{x~D}  (e^{-y F_{t-1} (x)} (yf_t(x) )) \\ &= \underset {f_t}{argmax} E_{x~D}  (\frac{ e^{-y F_{t-1} (x)}}{ E_{x~D} ( e^{-y F_{t-1}(x)}) } (yf_t(x) )) \\ &=\underset {f_t}{argmax} E_{x~D_t}  (yf_t(x))\\ &=\underset {f_t}{argmax} P(y=f_t(x))yf_t(x)+P(y \neq f_t(x))yf_t(x)\\ &=\underset {f_t}{argmax} P(y=f_t(x))yf_t(x)+(1- P(y=f_t(x)))yf_t(x)\\ &=\underset {f_t}{argmax} P(y=f_t(x))-(1- P(y=f_t(x)))\\ &=\underset {f_t}{argmax} 2P(y=f_t(x))-1 \end{align}

所以由这个公式可以看出,在用上述方法构造的D_t下,最大化t步分类器的正确率(最小化t步分类器误差),就是最小化最初的前t步的期望损失函数。所以这样构造的D_t是符合逻辑的。

D_t的通项表达式,可以得到D_t的递推表达公式:

\begin{align} D_{t+1}(x) &= \frac{ E_{x~D} ( e^{-y F_{t-1}(x)}) e^{-y F_{t} (x)}}{ E_{x~D} ( e^{-y F_{t}(x)}) e^{-y F_{t-1} (x)} }D_{t}(x)\\ &= \frac{ E_{x~D} ( e^{-y F_{t-1}(x)})}{ E_{x~D} ( e^{-y F_{t}(x)}) }D_{t}(x) e^{-y \alpha_t f_{t} (x)} \end{align}

这个公式中,\frac{ E_{x~D} ( e^{-y F_{t-1}(x)})}{ E_{x~D} ( e^{-y F_{t}(x)}) }是一个常数,对于同一个弱分类器,所有的样本x都一样,所有在更新权重时只需要计算D_{t}(x) e^{-y \alpha_t f_{t} (x)}然后再将其归一化即可。可以看到,在更新权重时,与F(x)没有任何关系,仅与f(x)相关。

至此,Adaboost的两个权重就推导出来了。

下面给出Adaboost的算法流程:

————————————————————————————————————————————————
输入:训练集D={(x_1,y_1),x_2,y_2),...(x_m,y_m)};
           基学习算法\xi;
           训练轮数T.
过程:
1:D_1(x)=1/m
2:for t=1,2,...T do
3: f_t=\xi(D,D_t)(基于样本分布D_t在样本数据D里面训练即学习基模型)
4: \epsilon_t=P_{x \in D_t}(f(x) \neq y)(计算第t步的错误率)
5: if \epsilon_t>0.5 then break;(若错误率大于0.5,直接结束)
6: \alpha_t=\frac{1}{2}ln(\frac{1-\epsilon_t}{\epsilon_t})
7: D_{t+1}(x)=\frac{D_t(x)e^{-\alpha_t f_t(x)}}{Z_t}(Z_t为归一化常数)
8:end for
输出:sign(\sum_{t=1}^T(\alpha_t f_t(x)))
————————————————————————————————————————————————

Adaboost有一个神奇的地方,是Aaboost不容易产生过拟合,由上面的分析,bagging是减小方差,boosting是减小偏差,但是Adaboost做为boosting的一种,在减小偏差的同时,也减小了方差,所以Aaboost不容易产生过拟合,是Adaboost的牛逼的地方。

至于Adaboost在理论上为什么不容易过拟合,周志华教授的《Ensemble Methods: Foundations and Algorithms》的第33~34页给出了具体的证明。不在此陈述。

下面给出Adaboost在分类问题上的一个例子。

#-*- conding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

from sklearn.ensemble import AdaBoostClassifier#adaboost引入方法
from sklearn.tree import DecisionTreeClassifier
from sklearn.datasets import make_gaussian_quantiles#造数据
## 设置属性防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False
## 创建数据
X1, y1 = make_gaussian_quantiles(cov=2.,
                                 n_samples=200, n_features=2,
                                 n_classes=2, random_state=1)#创建符合高斯分布的数据集
X2, y2 = make_gaussian_quantiles(mean=(3, 3), cov=1.5,
                                 n_samples=300, n_features=2,
                                 n_classes=2, random_state=1)

X = np.concatenate((X1, X2))
y = np.concatenate((y1, - y2 + 1))
plot_step = 0.02
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
                     np.arange(y_min, y_max, plot_step))
#构建adaboost模型
bdt = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1),
                         algorithm="SAMME.R",#可以不写
                         n_estimators=200)
#数据量大的时候,可以增加内部分类器的树深度,也可以不限制树深
#max_depth树深,数据量大的时候,一般范围在10——100之间
#数据量小的时候,一般可以设置树深度较小,或者n_estimators较小
#n_estimators 迭代次数或者最大弱分类器数:200次
#base_estimator:DecisionTreeClassifier 选择弱分类器,默认为CART树
#algorithm:SAMME 和SAMME.R 。运算规则,后者是优化算法,以概率调整权重,迭代速度快,
#需要能计算概率的分类器支持
#learning_rate:0<v<=1,默认为1,正则项 衰减指数
#loss:linear、‘square’exponential’。误差计算公式:一般用linear足够
bdt.fit(X, y)

#预测
Z = bdt.predict(np.c_[xx.ravel(), yy.ravel()])
#设置维度
Z = Z.reshape(xx.shape)
## 画图
plot_colors = "br"
class_names = "AB"

plt.figure(figsize=(10, 5), facecolor='w')
#局部子图
plt.subplot(121)
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)
for i, n, c in zip(range(2), class_names, plot_colors):
    idx = np.where(y == i)
    plt.scatter(X[idx, 0], X[idx, 1],
                c=c, cmap=plt.cm.Paired,
                label=u"class%s" % n)
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.legend(loc='upper right')
plt.xlabel('x')
plt.ylabel('y')
plt.title(u'AdaBoost accuracy:%.2f%%' % (bdt.score(X, y) * 100))

#获取决策函数的数值
twoclass_output = bdt.decision_function(X)
#获取范围
plot_range = (twoclass_output.min(), twoclass_output.max())
plt.subplot(122)
for i, n, c in zip(range(2), class_names, plot_colors):
#直方图
    plt.hist(twoclass_output[y == i],
             bins=20,
             range=plot_range,
             facecolor=c,
             label=u'class %s' % n,
             alpha=.5)
x1, x2, y1, y2 = plt.axis()
plt.axis((x1, x2, y1, y2 * 1.2))
plt.legend(loc='upper right')
plt.ylabel(u'Frequency')
plt.xlabel(u'F(x)')
plt.title(u'AdaBoost value')

plt.tight_layout()
plt.subplots_adjust(wspace=0.35)
plt.show()
/root/anaconda3/lib/python3.6/site-packages/matplotlib/font_manager.py:1297: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans
  (prop.get_family(), self.defaultFamily[fontext]))
output_180_1.png
AdaBoostClassifier与DecisionTreeRegressor调参

·max_features(划分时考虑的最大特征数): 可以使用很多种类型的值,默认是"None",意味着划分时考虑所有的特征数;如果是"log2"意味着划分时最多考虑log2N个特征;如果是"sqrt"或者"auto"意味着划分时最多考虑N−−√个特征。如果是整数,代表考虑的特征绝对数。如果是浮点数,代表考虑特征百分比,即考虑(百分比xN)取整后的特征数。其中N为样本总特征数。一般来说,如果样本特征数不多,比如小于50,我们用默认的"None"就可以了,如果特征数非常多,我们可以灵活使用刚才描述的其他取值来控制划分时考虑的最大特征数,以控制决策树的生成时间。

·max_depth(决策树最大深): 默认可以不输入,如果不输入的话,决策树在建立子树的时候不会限制子树的深度。一般来说,数据少或者特征少的时候可以不管这个值。如果模型样本量多,特征也多的情况下,推荐限制这个最大深度,具体的取值取决于数据的分布。常用的可以取值10-100之间。

·min_samples_split(内部节点再划分所需最小样本数): 这个值限制了子树继续划分的条件,如果某节点的样本数少于min_samples_split,则不会继续再尝试选择最优特征来进行划分。 默认是2.如果样本量不大,不需要管这个值。如果样本量数量级非常大,则推荐增大这个值。
· 叶子节点最少样本数min_samples_leaf: 这个值限制了叶子节点最少的样本数,如果某叶子节点数目小于样本数,则会和兄弟节点一起被剪枝。 默认是1,可以输入最少的样本数的整数,或者最少样本数占样本总数的百分比。如果样本量不大,不需要管这个值。如果样本量数量级非常大,则推荐增大这个值。

·min_weight_fraction_leaf(叶子节点最小的样本权重和):这个值限制了叶子节点所有样本权重和的最小值,如果小于这个值,则会和兄弟节点一起被剪枝。 默认是0,就是不考虑权重问题。一般来说,如果我们有较多样本有缺失值,或者分类树样本的分布类别偏差很大,就会引入样本权重,这时我们就要注意这个值了。

·max_leaf_nodes(最大叶子节点数): 通过限制最大叶子节点数,可以防止过拟合,默认是"None”,即不限制最大的叶子节点数。如果加了限制,算法会建立在最大叶子节点数内最优的决策树。如果特征不多,可以不考虑这个值,但是如果特征分成多的话,可以加以限制,具体的值可以通过交叉验证得到。

Adaboost应用

Adaboost可以用在人脸识别领域,具体算法推导和流程,在下一篇文章中介绍,在此不展开。

基于回归的Adaboost算法

对于Adaboost的回归算法,我未找到相关推导,只找到了算法流程,因为不明白其推导过程,就不在此陈述。

GBDT

GBDT的全名是Gradient Boosting Decision Tree,所以GBDT可以分成两部分学习,Gradient Boosting梯度提升,Decision Tree决策树,决策树已经在上文中进行讨论,一般使用CART(分类回归树),回归问题中,分裂时使最小二乘法,分类问题中,分裂是使用基尼系数。所以只要搞清楚Gradient Boosting就可以完全理解GBDT。

GBDT与Adaboost一样的地方都是一种提升算法,都是串行,下一步是根据上一步的结果来做判断,不一样的是,Adaboost是改变样本的权重,提高上一步分错样本的权重来训练下一步弱分类器,而GBDT是不改变样本,但是改变的是样本的标签y值,其第t步训练弱分类器时,目标变量(样本的标签)是上一步训练结束得到的残差,也就是对残差的拟合,这样,最后得到的每一步的值直接相加得到最终值,而不是像adaboost一样,是每一步的加权和。

GBDT即可以做回归问题,也可以做二分类问题,还可以做多分类问题。其推导过程是统一,这一点和Adaboost很不一样,Adaboost在分类算法推导的时候,必须确定其损失函数,并且必须要是指数损失函数,其他损失函数是推导不出来的,而GBDT不需要这样,GBDT的推导可以不用损失函数的表达式,直接用符号代替即可。

GBDT 算法推导

对于每一种模型和算法,其推导时,都可以从损失函数出发。GBDT也不例外,设:
F(x)= \sum_{t=1}^T( f_t(x))
为GBDT最后的预测值,若是回归问题其就是预测值,若是分类问题,还要将其映射的具体的类别。

那么L(y,F(x))为GBDT的期望损失函数。其中:y为样本的真实值,x为样本值。

由于F(x)是加法模型,所以和adaboost一样,使用贪心算法,追求当前最优。而和Adaboost不一样的是,adaboost损失函数优化过程中,未知的由两部分,一部分是每一个弱分类器的权重,另一部分是每一次训练弱分类器时样本的权重,而GBDT损失函数优化过程中,未知的只有一个,是当前弱分类器训练时样本的标签。

设在GBDT第t步分类时,样本的标签为y_t,则在第t步时,使用贪心算法,追求当前最优,损失函数为:
L(y, F_t(x))= L(y,\sum_{k=1}^{t-1}( f_k(x)) + y_t) \\ = L(y,F_{t-1}(x) + y_t)

上式由泰勒展开得到:
L(y,F_{t-1}(x) + y_t) = L(y,F_{t-1}(x)) + y_t \frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)}

由上式可以得到,当y_t=- \frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} 时,

L(y,F_{t-1}(x) + y_t) = L(y,F_{t-1}(x)) -(\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} )^2

所以当y_t=- \frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} 时,第t步后损失函数一定比第t-1步的损失函数要小,这样一直到T步结束,损失函数就会一直变小。

这个地方有个小问题,就是在回归问题里面第t步时,为什么不直接用真实值和前t-1步的预测值的残差y_t=y-F_{t-1}(x),而要费尽心思的去求负梯度去逼近残差,其实完全可以用残差,但是这里不使用直接使用残差,主要原因有两点,其一是,当损失函数取平方损失函数是,用负梯度代替残差的方式完全等价与直接用残差,因为直接用残差时,有\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)}=-(y-F_{t-1}(x)), 可以看到残差就是负的梯度方向。其二是,弱直接使用残差,相当于固定了损失函数为均方损失,不能用其他损失函数了,GBDT的牛比之处在于可以对任意一阶可导的损失函数套用,指数损失,log损失都可以 。

注:XGBoost的牛比之处在于可以对任意二阶可导的损失函数套用,指数损失,log损失函数都可以 。

由上面GBDT的推导可以看出,GBDT每一步是去拟合残差,也就是说每一步是去逼近一个值(残差),只有回归树才能去毕竟一个值,分类树是不能去逼近一个值的,所以,GBDT不管是处理分类问题还是回归问题,每一步使用的都是回归树。

下面给出GBDT的算法流程。

GBDT的算法流程

————————————————————————————————————————————————
输入:训练集D={(x_1,y_1),x_2,y_2),...(x_m,y_m)};
           基学习算法\xi;
           训练轮数T.
过程:
1:F_0(x)=argmin_c\sum_{i=1}^m L(y_i,c)初始化每一个样本对应的标签值
2:for t=1,2,...T do
3: \overline{y_i}=- \frac{\partial L(y_i, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x_i)}(计算每个点负梯度方向代替残差)
4: f_t=\xi(D,\overline{y})(基于样本数据D,与样本新的标签训练学习基模型,这里用到了叶子节点的最佳负梯度拟合的线性搜索)
5: F_t(x)=F_{t-1}+f_t(x)(迭代新的预测值)
6:end for
输出:F(x)=\sum_{t=1}^Tf_t(x)
————————————————————————————————————————————————

GBDT在处理回归问题,二分类问题,多分类问题时,虽然算法推导和算法流程是一致的,但在各个问题中,损失函数L(y,F(x))是不一致的,每个问题都有自己特定的期望损失函数L(y,F(x)).为什么可以使用这些各种各样的损失函数,本文不做介绍,会在专门的损失函数一文中,介绍这些损失函数的合理性和是怎么构造的。

GBDT回归问题

GBDT处理回归问题时,损失函数一般采用平方损失函数(MSE):
L(y,F(x))=\frac{1}{2}(y-F(x))^2

此时拟合的梯度如下:
\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} =F_{t-1}(x)-y

也可以用绝对误差(MAE):

L(y,F(x))=|y-F(x)|

其求导时,剔除不可导的点。此时拟合的梯度如下:
\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} =sign(y-F_{t-1}(x))

GBDT二分类问题

GBDT处理二分类问题时,损失函数一般采用对数似然损失函数:

L(y,F(x))=log(1+e^{-yF(x)})

此时拟合的梯度如下:
\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{t-1}(x)} =\frac{-y}{1+e^{yF_{t-1}(x)}}

此时,由于而二分类问题,使用的损失函数是对数似然损失函数,所以在算法的第4步,训练基模型去拟合负梯度时,与回归问题有点差别,因为回归问题的损失函数简单,求解基模型方便,而二分类用对数似然损失函数,导致在求解基模型时,优化问题不好解,一般采用近视的方法。具体如下:

在得到负梯度后,上述算法第4步,训练基模型:

c_{tj}=argmin_c\sum_{xi\in R_{tj}}L(y_i,F_{t-1}(x_i)+c) =argmin_c\sum_{xi\in R_{tj}}log(1+e^{-y_i(F_{t-1}(x_i)+c)})

其中c_{tj}为第t步第j个叶节点的值,R_{tj}是第t步第j个叶节点上的样本的集合。

此时,c_{tj}求解起来非常复杂。所以用一下公式代替。
c_{tj}=\frac{\sum_{x_i \in R_{tj} }r_{ti}}{\sum_{x_i \in R_{tj}} |r_{ti}|(1-|r_{ti}|)}

其中r_{ti}为第t步第i个样本点的负梯度值。既:
r_{ti} =\frac{y_i}{1+e^{y_iF_{t-1}(x_i)}}

GBDT多分类问题

GBDT 多分类问题用到的损失函数是softmax_cross_entropy,既:
L(y,F(x))=-\sum_{k=1}^Ky_klog(p_k(x))

其中:K是多分类问题的类别个数,y_k要么是0,要么是1,当x在第k_0个类别中时,y_{k_0}=1,其他的y_k=0.上式中的p_k(x),是softmax,表达式如下:
p_k(x)=\frac{e^{F_k(x)}}{\sum_{l=1}^Ke^{F_l(x)}}

其中:要想知道F_k(x)是什么,必须要理解,GBDT在K分类的时候,其实每一步是进行了K个二分类,在每一步训练K个基模型,F_k就是当前这一步训练的第k个模型的值。
所以这也可以看出GBDT实际训练了TxK个基模型。

所以上式还有个简单的表示方式:
L(y,F(x))=-log(p_k(x))

所以拟合的梯度如下:
\frac{\partial L(y, \chi)}{\partial\chi } |_{\chi=F_{l,t-1}(x)} =- y_{i,l}+p_{l,t-1}(x)

上式中y_{i,l}的意思是x_i属于l类别是1,不属于是零。

和二分类一样,用上述负梯度去训练下一步的模型,求模型叶节点上的值时,优化比较复杂,所以,第t步模型,第l个类别为1,其他类别为零,训练得到的决策树上第j个叶节点上的值为:
c_{t,j,l}=\frac{K-1}{K}\frac{\sum_{x_i \in R_{tjl} }r_{til}}{\sum_{x_i \in R_{tjl}} |r_{til}|(1-|r_{til}|)}

其中: r_{til}= y_{i,l}-p_{l,t-1}(x)

至此,GBDT的算法流程全部得到。

下面给出GBDT的一个例子。

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn import metrics

# 导入鸢尾花的数据
iris = datasets.load_iris()
# 特征数据
data = iris.data[:100] # 有4个特征
# 标签
label = iris.target[:100]

# 提取训练集和测试集
# random_state:是随机数的种子。
train_x, test_x, train_y, test_y = train_test_split(data, label,test_size=0.3, random_state=0)

# 模型训练,使用GBDT算法
gbr = GradientBoostingClassifier(n_estimators=3000, max_depth=2, min_samples_split=2, learning_rate=0.1)
gbr.fit(train_x, train_y)

y_pred = gbr.predict(test_x)

# ROC曲线下与坐标轴围成的面积
print ('AUC: %.4f' % metrics.roc_auc_score(test_y,y_pred))
# 准确率
print ('ACC: %.4f' % metrics.accuracy_score(test_y,y_pred))
print ('Recall: %.4f' % metrics.recall_score(test_y,y_pred))
# 精确率和召回率的调和平均数
print ('F1-score: %.4f' %metrics.f1_score(test_y,y_pred))
print ('Precesion: %.4f' %metrics.precision_score(test_y,y_pred))
metrics.confusion_matrix(test_y,y_pred)
AUC: 1.0000
ACC: 1.0000
Recall: 1.0000
F1-score: 1.0000
Precesion: 1.0000





array([[15,  0],
       [ 0, 15]])
GBDT参数

·n_estimators: 也就是弱学习器的最大迭代次数,或者说最大的弱学习器的个数。一般来说n_estimators太小,容易欠拟合,n_estimators太大,又容易过拟合,一般选择一个适中的数值。默认是100。在实际调参的过程中,我们常常将n_estimators和下面介绍的参数learning_rate一起考虑。

·learning_rate: 即每个弱学习器的权重缩减系数νν,也称作步长,在原理篇的正则化章节我们也讲到了,加上了正则化项,我们的强学习器的迭代公式为fk(x)=fk−1(x)+νhk(x)fk(x)=fk−1(x)+νhk(x)。νν的取值范围为0<ν≤10<ν≤1。对于同样的训练集拟合效果,较小的νν意味着我们需要更多的弱学习器的迭代次数。通常我们用步长和迭代最大次数一起来决定算法的拟合效果。所以这两个参数n_estimators和learning_rate要一起调参。一般来说,可以从一个小一点的νν开始调参,默认是1。

·subsample: 即我们在原理篇的正则化章节讲到的子采样,取值为(0,1]。注意这里的子采样和随机森林不一样,随机森林使用的是放回抽样,而这里是不放回抽样。如果取值为1,则全部样本都使用,等于没有使用子采样。如果取值小于1,则只有一部分样本会去做GBDT的决策树拟合。选择小于1的比例可以减少方差,即防止过拟合,但是会增加样本拟合的偏差,因此取值不能太低。推荐在[0.5, 0.8]之间,默认是1.0,即不使用子采样。

·init: 即我们的初始化的时候的弱学习器,拟合对应原理篇里面的f0(x)f0(x),如果不输入,则用训练集样本来做样本集的初始化分类回归预测。否则用init参数提供的学习器做初始化分类回归预测。一般用在我们对数据有先验知识,或者之前做过一些拟合的时候,如果没有的话就不用管这个参数了。

·loss: 即我们GBDT算法中的损失函数。分类模型和回归模型的损失函数是不一样的。

对于分类模型,有对数似然损失函数”deviance”和指数损失函数”exponential”两者输入选择。默认是对数似然损失函数”deviance”。在原理篇中对这些分类损失函数有详细的介绍。一般来说,推荐使用默认的”deviance”。它对二元分离和多元分类各自都有比较好的优化。而指数损失函数等于把我们带到了Adaboost算法。
对于回归模型,有均方差”ls”, 绝对损失”lad”, Huber损失”huber”和分位数损失“quantile”。默认是均方差”ls”。一般来说,如果数据的噪音点不多,用默认的均方差”ls”比较好。如果是噪音点较多,则推荐用抗噪音的损失函数”huber”。而如果我们需要对训练集进行分段预测的时候,则采用“quantile”。

·alpha:这个参数只有GradientBoostingRegressor有,当我们使用Huber损失”huber”和分位数损失“quantile”时,需要指定分位数的值。默认是0.9,如果噪音点较多,可以适当降低这个分位数的值。

XGBoost

正在拜读陈天奇大神的paper,发现里面细节非常多,打算真正搞懂以后再整理这一部分。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,117评论 4 360
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,963评论 1 290
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 107,897评论 0 240
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,805评论 0 203
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,208评论 3 286
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,535评论 1 216
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,797评论 2 311
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,493评论 0 197
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,215评论 1 241
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,477评论 2 244
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,988评论 1 258
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,325评论 2 252
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,971评论 3 235
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,055评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,807评论 0 194
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,544评论 2 271
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,455评论 2 266

推荐阅读更多精彩内容