ML算法：XGBoost
XGBoost
实验目的
本实验要求学习XGBoost模型的原理，并实现该模型，以及在数据集上完成训练和测试。
实验原理
XGBoost的基本原理就是以决策树作为基学习器的加法模型，每棵树根据前t个模型的输出计算节点的权重，最终的输出为：
本实验的问题模型为回归问题，即$𝑙𝑜𝑠𝑠(𝑦_𝑖,𝑦 _𝑖^{(𝑡)} )=(𝑦_𝑖−\hat y _𝑖^{(𝑡)} )^2$。
对于损失loss，记 $g_i$ 为一阶导数，$h_i$ 为二阶导数，即$g_i=\frac{\partial\ loss(y_i,y _i^{(t-1)})}{\partial\ y _i^{(t-1) } }=2(\hat y^{t-1}-y)$, $h_i=\frac{\partial^2 loss(y_i,y _i^{(t-1)} )}{\partial \ (y _i^{(t-1)} )^2 }=2\\ $。
在训练第t棵决策树时,目标函数为$Obj^{(t)}=∑_{i=1}^n[g_i f_t (x_i )+\frac12 h_i f_t^2 (x_i)]+penalty(f_t )$
对于一棵决策树，其惩罚为：
将分配到第 $j$ 个叶子节点的样本用 $I_j$ 表示，即 $I_j=\{i|q(x_i )=j\} (1≤j≤T)$，在树结构确定时可以进行如下优化：
求解$w$ :$w_j^*=-\frac{G_j}{H_j+\lambda}$，此时决策树的得分：$Obj^{(t)}=-\frac12\sum^{T}_{j=1}\frac{G_j^2}{H_j+\lambda}+\gamma T$
实验步骤
读取数据和数据处理
数据文件中的各项数据均不含有NULL项，无需进行特殊处理。
df = pd.read_csv('train.data')
df.shape[1]
df = pd.read_csv('train.data',names=[i for i in range(df.shape[1])])
#重新读取数据，防止将第一行的数据识别为索引
Data=np.array(df)
a=int(0.3*Data.shape[0])
Data=np.c_[Data,np.zeros(Data.shape[0]).T]
num=np.random.permutation(Data.shape[0])
Data=Data[num]
test_data=Data[0:a-1,:]
train_data=Data[a:,:]
#按照7：3的比例划分训练集和测试集在数据集的基础上添加一列0，用来存储每轮训练结束后的预测值，即data=[:,-2]代表真实值,data[:,-1]代表预测值。
回归树模型
超参数:
maxdepth=4,minscore=0,minindex=3,lambda=800,gamma=1e-7
class RegressionTree:
def __init__(self) :
...初始化
def Leaf(self,data):#叶节点
G=(2*np.sum(data[:,-1]-data[:,-2]))
H=2*data.shape[0]
return -G/(H+self.lam)
def Obj(self,data):#计算得分
G=(2*np.sum(data[:,-1]-data[:,-2]))
H=2*data.shape[0]
obj=-G**2*0.5/(H+self.lam)+self.gamma
return obj
def split(self,data,num,val):#根据num和val找到最佳划分
lchild=data[np.nonzero(data[:,num]<=val)[0],:]
rchild=data[np.nonzero(data[:,num]>val)[0],:]
return lchild,rchild
def ChooseSplit(self,data,depth):
"""
对样本集进行划分，找到最佳划分
"""
if data.shape[0]<self.minindex or depth>self.maxdepth:
return -1,self.Leaf(data)
if len(set(data[:,-2].T.tolist()))==1:
return -1,self.Leaf(data)
pobj=self.Obj(data)
maxscore=-1
maxnum=0
maxval=0
n,m=data.shape
for i in range(m-2):
n=data[:,i].T.tolist()
for j in set(n):
lchild,rchild=self.split(data,i,j)
if len(lchild)==0 or len(rchild)==0:
continue
score=-self.Obj(lchild)+self.Obj(rchild)+pobj
if score>maxscore:
maxscore=score
maxnum=i
maxval=j
if maxscore<=self.minscore:
return -1,self.Leaf(data)
return maxnum,maxval
def CreateTree(self,data,depth):
num,val=self.ChooseSplit(data,depth)
if num==-1:
return val
lchild,rchild=self.split(data,num,val)
Tree={}
"""
如果不是叶节点，则以字典形式建树，否则节点只存储叶子的值
"""
Tree['num']=num
Tree['val']=val
Tree['lchild']=self.CreateTree(lchild,depth+1)
Tree['rchild']=self.CreateTree(rchild,depth+1)
return Tree其中，决策树的非叶子节点有属性：
- num，val：在该结点处选取第num个属性值进行划分，最佳划分点为val
- lchild，rchild：分别指向第num个属性小于等于val和大于val的子树
对于单棵决策树的停止策略：
- maxdepth：每次划分记录树的高度，深度为maxdepth时直接当作叶子节点
- minindex：每次划分时，如果当前节点的样本数小于minindex，则直接作为叶子节点
- minscore：若每次划分的最佳划分的得分小于等于minscore时，直接作为叶子节点
优化ChooseSpilt函数。为了减少每次对得分Obj的计算时存在的重复计算，将函数的核心部分做以下修改：（将数据按照当前属性进行排序，每次选取的属性值依次增大，因此左右两孩子的g，h线性增/减，对于左孩子：$\Delta g=\sum _{data[j,i]=val}(data1[j,-1]-data1[j,-2])2.0$,$\Delta h=2\sum_{data[j,i]=val}$）
gl=0
hl=0
gr=(2*np.sum(data[:,-1]-data[:,-2]))
hr=2*n
data1=data[data[:,i].argsort()]
for j in range(n-1):
gl+=(data1[j,-1]-data1[j,-2])*2.0
gr-=(data1[j,-1]-data1[j,-2])*2.0
hl+=2
hr-=2
if data1[j,i]==data1[j+1,i]:
continue
obj1=-gl**2*0.5/(hl+self.lam)+self.gamma
obj2=-gr**2*0.5/(hr+self.lam)+self.gamma
score=pobj-obj1-obj2
if score>maxscore:
maxscore=score
maxnum=i
maxval=data1[j,i]XGBoost模型
超参数：treenum=30
class XGBoost:
def __init__(self,data,treenum=30):
...
def isTree(self,tree):#判断节点的类型，非叶子节点为字典类型，叶子节点为float类型
return type(tree).__name__=='dict'
def find(self,data,i,tree):#对于单个样本，找到所在决策树的位置
if self.isTree(tree):
if data[i,tree["num"]]<=tree['val']:
return self.find(data,i,tree['lchild'])
else:
return self.find(data,i,tree['rchild'])
else:
return tree
def fit(self):
"""
训练函数，其中共训练了treenum棵树，调用RegressionTree，并计算损失和相关系数
"""
loss=[]
r2=[]
va=(np.var(self.data[:,-2]))
for i in range(self.treenum):
print(i)
tree=Rtree.CreateTree(self.data,1)
for j in range(self.length):
self.data[j,-1]+=self.find(self.data,j,tree)
r=np.sqrt(np.sum((self.data[:,-2]-self.data[:,-1])**2/self.length))
r2.append(1-r**2/va)
lo=np.sum((self.data[:,-2]-self.data[:,-1])**2)
loss.append(lo)
self.forest.append(tree)
return loss,r2
def predict(self,data):
m,n=data.shape
for i in range(self.treenum):
for j in range(m):
test_data[j,-1]+=self.find(data,j,self.forest[i])训练模型并实现损失函数可视化
Rtree=RegressionTree()
boost=XGBoost(train_data)
loss,r2=boost.fit()
boost.predict(test_data)
x = np.arange(1,len(loss)+1)
plt.plot(x,loss)
plt.xlabel(u"times")
plt.ylabel(u"loss-value")
plt.show()模型预测并评估结果
结果评估指标：
- $RMSE=\sqrt{\frac1n\sum_{i=1}^n{(y-\hat y)^2}}$
- $R2=1-\frac{RMSE}{var(y)}$
- 误差分布直方图
r=np.sqrt(np.sum((test_data[:,-2]-test_data[:,-1])**2/test_data.shape[0]))
r2=(1-r**2/np.var(test_data[:,-2]))
res=test_data[:,-1]-test_data[:,-2]
res=np.sort(res)
plt.hist(res, bins=30, rwidth=0.9, density=True)
plt.xlabel(u"data")
plt.ylabel(u"num")
plt.show()
print(test_data[:,-1]-test_data[:,-2])
print(r,r2)实验结果与分析
（调参过程略）
