参考了统计学习方法,西瓜书,Machine Learnig with python做的总结,所以不能作为教程,还包含自己用sklearn做的一些对比实验,原文是写在jupyter上的,这里是直接转为.md导过来的,所以格式有些问题,有些东西还待完善…
注意几点:连续特征处理,预测问题或者说回归问题(连续性目标特征)
决策树(Decision tree)
熵
熵表示随机变量不确定性的度量。离散随机变量X的概率分布为,P(X=xi)=pi,i=1,2,3...,nP(X=x_i)=p_i,i=1,2,3...,nP(X=xi)=pi,i=1,2,3...,n.则随机变量X的熵可以定义为:H(p)=−∑i=1npilog2(pi)H(p)=-\sum_{i=1}^np_ilog_2(p_i)H(p)=−i=1∑npilog2(pi)0≤H(p)≤log2(n)0\leq H(p)\leq log_2(n)0≤H(p)≤log2(n)熵越大,随机变量的不确定性就越大.当随机变量的取任何值概率都相等时,也就是pi=1np_i=\frac{1}{n}pi=n1时,熵最大.此时可以知道H(p)=−n×1n×log2(1n)=log2(n)H(p)=-n\times\frac{1}{n}\times log_2(\frac{1}{n})=log_2(n)H(p)=−n×n1×log2(n1)=log2(n).
条件熵
我们知道条件概率为P(Y=yi∣X=xi)P(Y=y_i|X=x_i)P(Y=yi∣X=xi),表示在已知X=xiX=x_iX=xi条件下Y=yiY=y_iY=yi的概率.\quad则条件熵定义为在已知随机变量XXX的条件下,随机变量Y的不确定性,表示为H(Y∣X)H(Y|X)H(Y∣X):H(Y∣X)=−∑i=1npiH(Y∣X=xi)H(Y|X)=-\sum_{i=1}^n p_iH(Y|X=x_i)H(Y∣X)=−i=1∑npiH(Y∣X=xi)
信息增益
特征A对于训练数据D的信息增益表示为g(D,A)g(D,A)g(D,A):g(D,A)=H(D)−H(D∣A)g(D,A)=H(D)-H(D|A)g(D,A)=H(D)−H(D∣A)H(D)表示原始数据分类的不确定性,H(D∣A)H(D|A)H(D∣A)表示特征A给定条件下数据集D分类的不确定性.g(D,A)g(D,A)g(D,A)就表示给定特征A后数据集D不确定性减小的程度.对于数据集DDD和特征AAA有:H(D)=−∑i=1K∣Ck∣∣D∣log2∣Ck∣∣D∣H(D)=-\sum_{i=1}^K\frac{|C_k|}{|D|}log_2\frac{|C_k|}{|D|}H(D)=−i=1∑K∣D∣∣Ck∣log2∣D∣∣Ck∣,这里∣Ck∣|C_k|∣Ck∣表示类别CkC_kCk的数量.H(D∣A)=∑i=1n∣Di∣∣D∣H(Di)=∑i=1n∣Di∣∣D∣∑k=1K∣Dik∣∣Di∣log2(∣Dik∣∣Di∣)H(D|A)=\sum_{i=1}^n\frac{|D_i|}{|D|}H(D_i)=\sum_{i=1}^n\frac{|D_i|}{|D|}\sum_{k=1}^K\frac{|D_{ik}|}{|D_i|}log_2(\frac{|D_{ik}|}{|D_i|})H(D∣A)=i=1∑n∣D∣∣Di∣H(Di)=i=1∑n∣D∣∣Di∣k=1∑K∣Di∣∣Dik∣log2(∣Di∣∣Dik∣)这里nnn表示特征AAA的可取值数量.
信息增益比
定义:特征A对训练数据集D的信息增益比gR(D,A)g_R(D,A)gR(D,A)定义为其信息增益g(D,A)g(D,A)g(D,A)与训练数据D关于特征A的值的熵HA(D)H_A(D)HA(D)值比,即:gR(D,A)=g(D,A)HA(D)g_R(D,A)=\frac{g(D,A)}{H_A(D)}gR(D,A)=HA(D)g(D,A),其中HA(D)=−∑i=1n∣Di∣∣D∣log2∣Di∣∣D∣H_A(D)=-\sum_{i=1}^n\frac{|D_i|}{|D|}log_2\frac{|D_i|}{|D|}HA(D)=−∑i=1n∣D∣∣Di∣log2∣D∣∣Di∣,nnn是特征AAA取值个数.
基尼指数
定义:分类问题假设有KKK个类,样本点属于第kkk类的概率为pkp_kpk,则概率分布的基尼指数定义为Gini(p)=∑k=1Kpk(1−pk)=1−∑k=1Kpk2Gini(p)=\sum_{k=1}^Kp_k(1-p_k)=1-\sum_{k=1}^Kp^2_kGini(p)=k=1∑Kpk(1−pk)=1−k=1∑Kpk2
对于二分类问题,若样本点属于第一个类的概率是p,则概率分布的基尼指数为Gini(p)=2p(1−p)Gini(p)=2p(1-p)Gini(p)=2p(1−p)
对于给定样本集合DDD,其基尼指数为:Gini(D)=1−∑k=1K(CkD)2(1)Gini(D)=1-\sum_{k=1}^K(\frac{C_k}{D})^2 (1)Gini(D)=1−k=1∑K(DCk)2(1)
这里,CkC_kCk是DDD中属于第kkk类的样本自己,KKK是类的个数.
如果样本集合根据特征AAA是否取某一可能值aaa,被分割成D1D_1D1和D2D_2D2两个部分,即D1={(x,y)∈D∣A(x)=a},D2=D−D1D_1=\{(x,y)\in D|A(x)=a\},D_2=D-D_1D1={(x,y)∈D∣A(x)=a},D2=D−D1
则在特征AAA的条件下,集合DDD的基尼指数定义为:Gini(D,A)=D1DGini(D1)+D2DGini(D2)Gini(D,A)=\frac{D_1}{D}Gini(D_1)+\frac{D_2}{D}Gini(D_2)Gini(D,A)=DD1Gini(D1)+DD2Gini(D2)
基尼指数Gini(D)Gini(D)Gini(D)表示集合DDD的不确定性,基尼指数G(D,A)G(D,A)G(D,A)表示经A=aA=aA=a分割后集合DDD的不确定性.基尼指数越大,样本集合的不确定性也就越大.怎么理解:当种类数变多时(1)式中第二项会变小,这样会使基尼指数变大,这样也就是表示了不确定性越大.
决策树生成
ID3:
输入:训练数据集DDD,特征集AAA,阈值ϵ\epsilonϵ;
输出:决策树TTT
1.若DDD中所有实例属于同一类CkC_kCk,则TTT为单节点树,并将类CkC_kCk作为该节点的类标记,返回T;
2.若A=∅A=\varnothingA=∅,则返回单节点树,并将DDD中实例数最大的类CkC_kCk作为该节点的类标记,返回TTT;
3.否则,选择信息增益最大的特征AgA_gAg;
4.如果AgA_gAg信息增益小于阈值ϵ\epsilonϵ,则置TTT为单节点树,并将DDD中实例数最大的类CkC_kCk作为该节点的类标记,返回TTT;
5.否则,对AgA_gAg的每个可能值aia_iai,依Ag=aiA_g=a_iAg=ai将DDD分为若干个非空子集DiD_iDi,将DiD_iDi中实例数最大的类作为标记,构建子节点,由节点及其子节点构成树TTT,返回TTT;
6.对第i个子节点,以DiD_iDi为训练集,以A−{Ag}A-\{A_g\}A−{Ag}为特征集,递归调用1-5步,得到子树,返回TiT_iTi.
C45:
C45与ID3算法相似,但是使用的是信息增益比来选择特征,这是因为 以信息增益比作为划分训练数据集的特征,存在偏向于选择取值较多的特征的问题.使用信息增益比可以对这一问题进行校正. 这是因为AAA的取值很多的话∣Di∣∣D∣\frac{|D_i|}{|D|}∣D∣∣Di∣会很小,这样会导致H(D∣A)H(D|A)H(D∣A)也会很小,这样信息增益就会很大.而信息增益比在信息增益的基础上除以了特征AAA的取值类别的熵,也就是HA(D)H_A(D)HA(D)就表示特征AAA的取值不确定性程度或者说取值复杂度,那么信息增益的基础上除以HA(D)H_A(D)HA(D)就可以抵消掉特征AAA的取值复杂度的影响.
剪枝
决策树剪枝通过最小化决策树整体的损失函数来实现.
CART
后面回归树一起
优点
不需要特征归一化处理;
可以做分类,也可以做回归;
可以建模非线性关系;
缺点
若是连续特征,树可能会相当大;
决策树的一个缺点是要求训练数据中包含所有标签类,否则对于数据中没有出现过的标签类没有判别能力
数据微小变化可能导致不同的树;
若特征很多,而数据量少,那么树很容易过拟合;
python实现
import pandas as pd
import numpy as np
from pprint import pprint
dataset=pd.read_csv("data/zoo.data",names=["animal_name","hair","feathers","eggs","milk",
"airbone","aquatic","predator","toothed","backbone","breathes","venomous","fins","legs"
,"tail","domestic","catsize","class"])##若有数据没有columns可以使用names来添加
dataset.head()
animal_name | hair | feathers | eggs | milk | airbone | aquatic | predator | toothed | backbone | breathes | venomous | fins | legs | tail | domestic | catsize | class | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | aardvark | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 0 | 0 | 1 | 1 |
1 | antelope | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 4 | 1 | 0 | 1 | 1 |
2 | bass | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 4 |
3 | bear | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 0 | 0 | 1 | 1 |
4 | boar | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 1 | 0 | 1 | 1 |
dataset=dataset.drop("animal_name",axis=1)#去掉列
dataset.head()
#dataset.loc[set(range(0,101))-set([9,100,2])]#去掉行
#dataset[["hair","eggs"]].iloc[[1,2,3,6]]#去掉行和列
#np.unique(dataset[["hair","eggs"]].iloc[[1,2,3,6]])#去掉重复项
#p.unique(dataset["tail"])
hair | feathers | eggs | milk | airbone | aquatic | predator | toothed | backbone | breathes | venomous | fins | legs | tail | domestic | catsize | class | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 0 | 0 | 1 | 1 |
1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 4 | 1 | 0 | 1 | 1 |
2 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 4 |
3 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 0 | 0 | 1 | 1 |
4 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 4 | 1 | 0 | 1 | 1 |
def entropy(target_col):
elements,counts = np.unique(target_col,return_counts = True)
entropy = np.sum([(-counts[i]/np.sum(counts))*np.log2(counts[i]/np.sum(counts)) for i in range(len(elements))])
return entropy
def InfoGain(data,split_attribute_name,target_name="class"):
total_entropy = entropy(data[target_name])
vals,counts= np.unique(data[split_attribute_name],return_counts=True)
Weighted_Entropy = np.sum([(counts[i]/np.sum(counts))*entropy(data.where(data[split_attribute_name]==vals[i]).dropna()[target_name]) for i in range(len(vals))])
Information_Gain = total_entropy - Weighted_Entropy
return Information_Gain
def ID3(data,originaldata,features,target_attribute_name="class",parent_node_class = None):
if len(np.unique(data[target_attribute_name])) <= 1:#也就是只有一个类别
return np.unique(data[target_attribute_name])[0]#返回该类
elif len(data)==0:
return np.unique(originaldata[target_attribute_name])[np.argmax(np.unique(originaldata[target_attribute_name],return_counts=True)[1])]
elif len(features) ==0:
return parent_node_class
else:
parent_node_class = np.unique(data[target_attribute_name])[np.argmax(np.unique(data[target_attribute_name],return_counts=True)[1])]
item_values = [InfoGain(data,feature,target_attribute_name) for feature in features] #Return the information gain values for the features in the dataset
best_feature_index = np.argmax(item_values)
best_feature = features[best_feature_index]
tree = {best_feature:{}}
features = [i for i in features if i != best_feature]
for value in np.unique(data[best_feature]):
value = value
sub_data = data.where(data[best_feature] == value).dropna()
subtree = ID3(sub_data,dataset,features,target_attribute_name,parent_node_class)
tree[best_feature][value] = subtree
#print(tree)
return(tree)
def predict(query,tree,default = 1):
for key in list(query.keys()):
if key in list(tree.keys()):
#2.
try:
result = tree[key][query[key]]
except:
return default
#3.
result = tree[key][query[key]]
#4.
if isinstance(result,dict):
return predict(query,result)
else:
return result
def train_test_split(dataset):
training_data = dataset.iloc[:80].reset_index(drop=True)
testing_data = dataset.iloc[80:].reset_index(drop=True)
return training_data,testing_data
training_data = train_test_split(dataset)[0]
testing_data = train_test_split(dataset)[1]
tree = ID3(training_data,training_data,training_data.columns[:-1])
pprint(tree)
{'legs': {0: {'fins': {0.0: {'toothed': {0.0: 7.0, 1.0: 3.0}},
1.0: {'eggs': {0.0: 1.0, 1.0: 4.0}}}},
2: {'hair': {0.0: 2.0, 1.0: 1.0}},
4: {'hair': {0.0: {'toothed': {0.0: 7.0, 1.0: 5.0}}, 1.0: 1.0}},
6: {'aquatic': {0.0: 6.0, 1.0: 7.0}},
8: 7.0}}
def test(data,tree):
#Create new query instances by simply removing the target feature column from the original dataset and
#convert it to a dictionary
queries = data.to_dict(orient = "records")
#Create a empty DataFrame in whose columns the prediction of the tree are stored
predicted = pd.DataFrame(columns=["predicted"])
#Calculate the prediction accuracy
for i in range(len(data)):
predicted.loc[i,"predicted"] = predict(queries[i],tree,1.0)
print('The prediction accuracy is: ',(np.sum(predicted["predicted"] == data["class"])/len(data))*100,'%')
testing_data.to_dict(orient = "records")#orient若不设置,默认key值是数值
test(testing_data,tree)
The prediction accuracy is: 85.71428571428571 %
sklearn实现
from sklearn.tree import DecisionTreeClassifier
dataset1=pd.read_csv("data\zoo.data",names=["animal_name","hair","feathers","eggs","milk",
"airbone","aquatic","predator","toothed","backbone","breathes","venomous","fins","legs"
,"tail","domestic","catsize","class"])
dataset1=dataset1.drop("animal_name",axis=1)
train_features=dataset1.iloc[:80,:-1]
train_features
test_features=dataset1.iloc[80:,:-1]
train_targets=dataset1.iloc[:80,-1]
test_targets=dataset1.iloc[80:,-1]
test_targets
80 3
81 7
82 4
83 2
84 1
85 7
86 4
87 2
88 6
89 5
90 3
91 3
92 4
93 1
94 1
95 2
96 1
97 6
98 1
99 7
100 2
Name: class, dtype: int64
model=DecisionTreeClassifier(criterion="entropy",max_depth=6,min_samples_leaf=1,min_samples_split=4)#gini
model.fit(train_features,train_targets)
#Scikit-Learn 用的是 CART 算法, CART 算法仅产生二叉树:每一个非叶节点总是只有
#两个子节点(只有是或否两个结果) 。然而,像 ID3 这样的算法可以产生超过两个子节,因此这里用sklearn没有python写的效果好
#点的决策树模型。
DecisionTreeClassifier(class_weight=None, criterion='entropy', max_depth=6,
max_features=None, max_leaf_nodes=None,
min_impurity_decrease=0.0, min_impurity_split=None,
min_samples_leaf=1, min_samples_split=4,
min_weight_fraction_leaf=0.0, presort=False, random_state=None,
splitter='best')
prediction=model.predict(test_features)
print("The prediction accuracy is:",model.score(test_features,test_targets)*100,"%")
The prediction accuracy is: 80.95238095238095 %
sklearn 中决策树的剪枝通过调参实现
clf = tree.DecisionTreeClassifier()这个构建决策树的构造函数,带有参数常用的包括如下:
criterion='gini', 选用基尼系数作为选择特征的分裂点“entropy”
max_depth=None, 树的最大深度
min_samples_split=2, 分裂点的样本个数
min_samples_leaf =1, 叶子节点的样本个数
max_leaf_nodes=None,最大的叶子节点数
### 试试随机森林
from sklearn.ensemble import RandomForestClassifier
model=RandomForestClassifier(n_estimators=5,max_leaf_nodes=6,n_jobs=-1,criterion="entropy")
model.fit(train_features,train_targets)
prediction=model.predict(test_features)
print("The prediction accuracy is:",model.score(test_features,test_targets)*100,"%")
The prediction accuracy is: 76.19047619047619 %
### 试试AdaBoostClassifier(或者说Booting DTC,提升决策分类树)
from sklearn.ensemble import AdaBoostClassifier
model=AdaBoostClassifier(DecisionTreeClassifier(max_depth=3),n_estimators=200,algorithm="SAMME.R",learning_rate=0.5)
model.fit(train_features,train_targets)
#model=AdaBoostClassifier(n_estimators=200,algorithm="SAMME.R",learning_rate=0.5)#SAMME.R
#model.fit(train_features,train_targets)
prediction=model.predict(test_features)
print("The prediction accuracy is:",model.score(test_features,test_targets)*100,"%")
The prediction accuracy is: 85.71428571428571 %
### 试试Gradient Tree Boosting(梯度提升决策分类树也就是常说的GBDT())
from sklearn.ensemble import GradientBoostingClassifier
model=GradientBoostingClassifier(max_depth=6,n_estimators=100,learning_rate=1)
model.fit(train_features,train_targets)
prediction=model.predict(test_features)
print("The prediction accuracy is:",model.score(test_features,test_targets)*100,"%")
The prediction accuracy is: 85.71428571428571 %
连续属性值的处理
对于连续属性值取值数目不在有限,因此不能直接根据连续属性的可取值进行划分。最简单的是采用二分法对连续属性进行处理,这正是C4.5决策树算法中采用的机制.见西瓜书P83
回归树(Regression Tree)
如果我们需要用树结构来做预测问题,比如使用属性房间数,地理位置来预测目标特征(target feature)即房屋价格,此时价格就是连续的。我们就需要使用回归树来解决这个问题。
回归树的生成与决策树生成基本一样,只是有两点改变,首先我们回顾一下决策树生成叶子节点时的停止准侧(标准critera):
1.如果拆分过程导致数据集为空,则返回原始数据的目标特征值
2.如果拆分过程使得数据无特征剩余,则返回父节点的目标特征值
3.如果拆分过程使得数据目标特征值是一致时,返回这个值
1.现在我们来考虑连续问题的情况,此时停止准则中的第三个点就不在适用,因为目标特征值是连续的,所以不可能拆分到一个纯的目标特征值.为了解决这个问题,我们可以使用一种个提前结束准则,即返回目标特征值的平均值,当拆分到数据集中数量小于等于5时.也就是在回归树中,我们采用平均目标特征值作为叶子结点(预测值)。注意这个5是可调的下面实验中我们将展示其影响
2.现在来考虑划分标准,我们希望通过这个划分标准划分得到的预测值,尽量靠近真实值,因此我们选取加权方差(Varience)作为划分标准。为甚不用熵,因为目标特征值很多情况下就一个,那么条件熵是0.
举个例子
WeightVar(Season)=19×(79−79)2+59×(352−211.8)2+(421−211.8)2+(12−211.8)2+(162−211.8)2+(112−211.8)24+19×(161−161)2+29×(109−137)2+(165−137)21=16429.1WeightVar(Season)=\frac{1}{9}\times(79-79)^2+\frac{5}{9}\times\frac{(352-211.8)^2+(421-211.8)^2+(12-211.8)^2+(162-211.8)^2+(112-211.8)^2}{4}+\frac{1}{9}\times(161-161)^2+\frac{2}{9}\times\frac{(109-137)^2+(165-137)^2}{1}=16429.1WeightVar(Season)=91×(79−79)2+95×4(352−211.8)2+(421−211.8)2+(12−211.8)2+(162−211.8)2+(112−211.8)2+91×(161−161)2+92×1(109−137)2+(165−137)2=16429.1
WeightVar(Weekday)=29×(109−94)2+(79−94)21+29×(162−137)2+(112−137)21+19×(421−421)2+29×(161−86.5)2+(12−86.5)21+29×(352−258.5)2+(165−258.5)21=6730WeightVar(Weekday)=\frac{2}{9}\times\frac{(109-94)^2+(79-94)^2}{1}+\frac{2}{9}\times\frac{(162-137)^2+(112-137)^2}{1}+\frac{1}{9}\times(421-421)^2+\frac{2}{9}\times\frac{(161-86.5)^2+(12-86.5)^2}{1}+\frac{2}{9}\times\frac{(352-258.5)^2+(165-258.5)^2}{1}=6730WeightVar(Weekday)=92×1(109−94)2+(79−94)2+92×1(162−137)2+(112−137)2+91×(421−421)2+92×1(161−86.5)2+(12−86.5)2+92×1(352−258.5)2+(165−258.5)2=6730
WeightVar(Weathersit)=49×(421−174.2)2+(165−174.2)2+(12−174.2)2+(161−174.2)2+(112−174.2)24+29×(352−230.5)2+(109−230.5)21+29×(79−120.5)2+(112−120.5)21=19646.83WeightVar(Weathersit)=\frac{4}{9}\times\frac{(421-174.2)^2+(165-174.2)^2+(12-174.2)^2+(161-174.2)^2+(112-174.2)^2}{4}+\frac{2}{9}\times\frac{(352-230.5)^2+(109-230.5)^2}{1}+\frac{2}{9}\times\frac{(79-120.5)^2+(112-120.5)^2}{1}=19646.83WeightVar(Weathersit)=94×4(421−174.2)2+(165−174.2)2+(12−174.2)2+(161−174.2)2+(112−174.2)2+92×1(352−230.5)2+(109−230.5)2+92×1(79−120.5)2+(112−120.5)2=19646.83
由于WeekdayWeekdayWeekday有最低的加权方差,因此选择这个特征作为根节点.
最后准确度度量使用的是根均方误差RMSE:RMSE=∑i=12(ti−Model(testi))2nRMSE=\sqrt{\frac{\sum_{i=1}^2(t_i-Model(test_i))^2}{n}}RMSE=n∑i=12(ti−Model(testi))2
CART
CART对分类树使用基尼指数最小化准则,进行特征选择,生成 二叉树.对于回归树采用平方误差最小化准则(没有使用RMSE) 由于CART生成二叉树,无论分类树还是回归树事实上对属性值选取采取的是二分法(连续属性处理方式).(ID3,C45不一定生成二叉树)
CART分类树
输入:训练数据集D,停止计算条件;
输出:CART决策树
(1) 设节点的训练数据集为DDD,对于每个特征AAA,对其可能的取的每个aaa值,计算A=aA=aA=a的基尼指数Gini(D,A=a).
(2) 在所有可能的特征以及它们所有可能的切分点aaa中,选择基尼指数最小的特征及其对应的切分点作为最优特征与最优切分点,生成两个子节点,并将数据集分配到两个子节点中。
(3) 对两个子节点递归调用(1),(2),直至满足停止条件
(4) 生成CART决策树
算法停止条件:节点中样本个数小于约定阈值,或样本集合的基尼指数小于预定阈值(样本基本属于同一类),或者没有更多特征.
CART回归树(最小二乘回归树)
CART回归树也叫最小二乘回归树,这是因为其以最小化平方误差来生成决策树的.其特征是采用二分法来选取的.而且其目标特征值实际上采用的就是平均处理的方式.
对于训练数据D={(x1,y1),(x2,y2),..,(xn,yn)}D=\{(x_1,y_1),(x_2,y_2),..,(x_n,y_n)\}D={(x1,y1),(x2,y2),..,(xn,yn)}
一个回归树对应着输入空间(特征空间)的一个划分,以及在划分单元上的的输出值.假设将输入空间划分为MMM个单元R1,R2,..,RMR_1,R_2,..,R_MR1,R2,..,RM,并且在每个单元RmR_mRm上有一个固定输出值cmc_mcm,于是回归树模型可以表示为f(x)=∑m=1McmI(x∈Rm)f(x)=\sum_{m=1}^Mc_mI(x\in R_m)f(x)=m=1∑McmI(x∈Rm)
那么这个回归树的误差可以表示∑i=1n(yi−ci)2\sum_{i=1}^n(y_i-c_i)^2i=1∑n(yi−ci)2
其中ci=f(x)∈{c1,c2,...,cM}c_i=f(x)\in\{c_1,c_2,...,c_M\}ci=f(x)∈{c1,c2,...,cM},由于回归树生成是最小平方误差,所以cmc_mcm应该是单元RmR_mRm中的所有输入实例xix_ixi对应的标签值yiy_iyi的均值,即c^m=cm=ave(yi∣xi∈Rm)\hat{c}_m=c_m=ave(y_i|x_i\in R_m)c^m=cm=ave(yi∣xi∈Rm)
这和上面讲的平均目标特征值作为叶子结点的处理方式是一样的
特征空间的划分采用启发式的方法,选择第jjj个变量(特征)和它的取值sss,作为切分变量和切分点,将数据集划分为两个.寻找最优切分变量和切分点通过最小化平方差.具体见算法
CART回归树(最小二乘回归树)算法
输入:训练数据集DDD;
输出:回归树f(x)f(x)f(x)
在训练数据集所在的输入空间中,递归地将每个区域划分为两个子区域,并决定每个子区域上的输出值,构建二叉树
(1)选择最优切分变量(特征)jjj与切分点sss,求解minj,s[minc1∑xi∈R1(j,s)(yi−c1)2+minc2∑xi∈R2(j,s)(yi−c2)2]min_{j,s}[min_{c_1}\sum_{x_i\in R_1(j,s)}(y_i-c_1)^2+min_{c_2}\sum_{x_i\in R_2(j,s)}(y_i-c_2)^2]minj,s[minc1xi∈R1(j,s)∑(yi−c1)2+minc2xi∈R2(j,s)∑(yi−c2)2]遍历变量jjj,对固定切分变量扫描切分点sss,选择使上式达到最小值得(j,s)(j,s)(j,s)
(2)用选定的对(j,s)划分区域,并决定相应的输出值:R1(j,s)={x∣xj≤s},R2(j,s)={x∣xj>s}R_1(j,s)=\{x|x^{j}\leq s\},R_2(j,s)=\{x|x^{j}>s\}R1(j,s)={x∣xj≤s},R2(j,s)={x∣xj>s}
c^m=1Nm∑xi∈Rm(j,s)yi,x∈Rm,m=1,2\hat{c}_m=\frac{1}{N_m}\sum_{x_i\in R_m(j,s)y_i},x\in R_m,m=1,2c^m=Nm1xi∈Rm(j,s)yi∑,x∈Rm,m=1,2
(3)继续对两个子区域调用(1),(2),直至满足停止条件.
(4)将输入空间划分为MMM个区域R1,R2,...,RMR_1,R_2,...,R_MR1,R2,...,RM,生成决策树:f(x)=∑m=1Mc^mI(x∈Rm)f(x)=\sum_{m=1}^M\hat{c}_mI(x\in R_m)f(x)=m=1∑Mc^mI(x∈Rm)
sklearn实现
from sklearn.tree import DecisionTreeRegressor
import pandas as pd
import numpy as np
dataset=pd.read_csv("data/Bike-Sharing-Dataset/day.csv")
dataset=dataset[['season','holiday','weekday','workingday','weathersit','cnt']].sample(frac=1)
#dataset.sample(frac=0.11)#随机抽取,参数frac:0-1之间,表示随机抽取比例,想抽取n个则直接设置参数n=2
#dataset.sample(n=5)
dataset.head()
season | holiday | weekday | workingday | weathersit | cnt | |
---|---|---|---|---|---|---|
706 | 4 | 0 | 5 | 1 | 2 | 5008 |
75 | 1 | 0 | 4 | 1 | 1 | 2744 |
263 | 3 | 0 | 3 | 1 | 2 | 4352 |
334 | 4 | 0 | 4 | 1 | 1 | 3727 |
416 | 1 | 0 | 2 | 1 | 1 | 3777 |
mean_data=np.mean(dataset.iloc[:,-1])#对最后一列取个平均
mean_data
4504.3488372093025
def train_test_split(dataset):
training_data=dataset.iloc[:int(0.7*len(dataset))].reset_index(drop=True)#drop index并重新设置
testing_data=dataset.iloc[int(0.7*len(dataset)):].reset_index(drop=True)
return training_data,testing_data
training_data=train_test_split(dataset)[0]
testing_data=train_test_split(dataset)[1]
#training_data
#testing_data.iloc[:,:-1]
regression_model=DecisionTreeRegressor(criterion="mse",min_samples_leaf=5)
regression_model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
max_leaf_nodes=None, min_impurity_decrease=0.0,
min_impurity_split=None, min_samples_leaf=5,
min_samples_split=2, min_weight_fraction_leaf=0.0,
presort=False, random_state=None, splitter='best')
predicted=regression_model.predict(testing_data.iloc[:,:-1])
predicted.shape
testing_data.iloc[:,-1:].shape
test_target=np.asarray(testing_data.iloc[:,-1:])
test_target=test_target.reshape([220,])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(test_target)-1))
RMSE
1579.8360862692232
现在我们看看:拆分数据集最小量变化对树结构的影响
import matplotlib.pyplot as pl
train_data_RMSE=[]
test_data_RMSE=[]
train_target=np.asarray(training_data.iloc[:,-1:]).reshape([len(training_data),])
test_taregt=np.asarray(testing_data.iloc[:,-1:]).reshape([len(testing_data,)])
for i in range(1,100):
model=DecisionTreeRegressor(criterion="mse",min_samples_leaf=i)
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
train_predict=model.predict(training_data.iloc[:,:-1])
test_predict=model.predict(testing_data.iloc[:,:-1])
train_rmse=np.sqrt(np.sum((train_predict-train_target)**2)/(len(training_data)-1))
test_rmse=np.sqrt(np.sum((test_predict-test_target)**2)/(len(testing_data-1)))
train_data_RMSE.append(train_rmse)
test_data_RMSE.append(test_rmse)
pl.plot(range(1,100),train_data_RMSE,label="train_data_RMSE")#也就是两条线相交的那个点是最好的
pl.plot(range(1,100),test_data_RMSE,label="test_data_RMSE")
pl.legend()
pl.show()
其他模型在此数据集上的表现:
from sklearn.linear_model import LinearRegression
model=LinearRegression(normalize=True)
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=True)
predicted=model.predict(testing_data.iloc[:,:-1])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(testing_data)-1))
print(RMSE)#发现方差更大,可能是因为数据并非是非线性的
33699.118262229465
#再试试Lasso和Ridge
from sklearn.linear_model import Lasso
model=Lasso(alpha=0.3,normalize=True)
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
Lasso(alpha=0.3, copy_X=True, fit_intercept=True, max_iter=1000,
normalize=True, positive=False, precompute=False, random_state=None,
selection='cyclic', tol=0.0001, warm_start=False)
predicted=model.predict(testing_data.iloc[:,:-1])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(testing_data)-1))
print(RMSE)#发现方差更大,可能是因为数据并非是非线性的
1710.6078966005025
#再试试Ridge
from sklearn.linear_model import Ridge
model=Ridge(0.3,normalize=True)
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
Ridge(alpha=0.3, copy_X=True, fit_intercept=True, max_iter=None,
normalize=True, random_state=None, solver='auto', tol=0.001)
predicted=model.predict(testing_data.iloc[:,:-1])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(testing_data)-1))
print(RMSE)#发现方差更大,可能是因为数据并非是非线性的
32156.18744278281
### 再试试AdaBoosting regression(提升回归树)
from sklearn.ensemble import AdaBoostRegressor
from sklearn.tree import DecisionTreeRegressor
model=AdaBoostRegressor(base_estimator=DecisionTreeRegressor(max_depth=3),n_estimators=100,learning_rate=1)#没有指出base_estimator,则默认是DecisionTreeRegressor(max_depth=3)
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
predicted=model.predict(testing_data.iloc[:,:-1])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(testing_data)-1))
print(RMSE)#发现方差更大,可能是因为数据并非是非线性的
1495.9294447250152
D:\anaconda\lib\site-packages\sklearn\utils\validation.py:761: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
y = column_or_1d(y, warn=True)
### 再试试Gradient Boosting Tree Regressor(梯度提升回归树)
from sklearn.ensemble import GradientBoostingRegressor
model=GradientBoostingRegressor(loss="huber",max_depth=2,n_estimators=200,learning_rate=0.05)#由于使用梯度,所以学习率更小,损失函数可选有ls,lad,huber,quamtile
model.fit(training_data.iloc[:,:-1],training_data.iloc[:,-1:])
predicted=model.predict(testing_data.iloc[:,:-1])
RMSE=np.sqrt(np.sum((predicted-test_target)**2)/(len(testing_data)-1))
print(RMSE)#发现方差更大,可能是因为数据并非是非线性的
D:\anaconda\lib\site-packages\sklearn\utils\validation.py:761: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
y = column_or_1d(y, warn=True)
1491.1037993981756