我正在尝试编写我自己的梯度提升算法。我知道有现成的软件包如 gbm 和 xgboost, 但我想通过自己编写来理解算法的工作原理。
我使用的是 iris 数据集，我的目标变量是 Sepal.Length（连续型）。我的损失函数是 mean(1/2*(y-yhat)^2)（基本上是均方误差前面乘以1/2），因此我的对应梯度只是残差 y - yhat。我将初始预测值设为0。
library(rpart)data(iris)#定义梯度函数grad.fun <- function(y, yhat) {return(y - yhat)}mod <- list()grad_boost <- function(data, learning.rate, M, grad.fun) { # 初始化拟合值为0 fit <- rep(0, nrow(data)) grad <- grad.fun(y = data$Sepal.Length, yhat = fit) # 初始化模型 mod[[1]] <- fit # 循环进行M次迭代 for(i in 1:M){ # 用梯度拟合基础学习器（树） tmp <- data$Sepal.Length data$Sepal.Length <- grad base_learner <- rpart(Sepal.Length ~ ., data = data, control = ("maxdepth = 2")) data$Sepal.Length <- tmp # 通过拟合当前模型得到拟合值 fit <- fit + learning.rate * as.vector(predict(base_learner, newdata = data)) # 更新梯度 grad <- grad.fun(y = data$Sepal.Length, yhat = fit) # 存储当前模型（索引是i + 1，因为i = 1包含了初始估计值） mod[[i + 1]] <- base_learner } return(mod)}有了这些，我将 iris 数据集分成训练和测试数据集，并将我的模型拟合到这些数据集上。
train.dat <- iris[1:100, ]test.dat <- iris[101:150, ]learning.rate <- 0.001M = 1000my.model <- grad_boost(data = train.dat, learning.rate = learning.rate, M = M, grad.fun = grad.fun)现在我计算 my.model 的预测值。对于 my.model，拟合值是 0（初始估计值向量）+ learning.rate * 树1的预测 + learning rate * 树2的预测 + ... + learning.rate * 树M的预测。
yhats.mymod <- apply(sapply(2:length(my.model), function(x) learning.rate * predict(my.model[[x]], newdata = test.dat)), 1, sum)# 计算RMSE> sqrt(mean((test.dat$Sepal.Length - yhats.mymod)^2))[1] 2.612972我有一些问题
- 我的梯度提升算法看起来正确吗？
- 我计算的预测值
yhats.mymod正确吗？
回答：
- 是的，看起来是正确的。在每个步骤中，你都在拟合伪残差，这些残差是损失相对于拟合值的导数。你在问题开始时正确地推导了这个梯度，甚至还正确处理了2的因子。
- 这也看起来是正确的。你按照训练时的方式，对模型进行加权聚合。
但对于未提及的问题，我想指出你的训练设置有几个奇怪之处。
iris数据集在3个种类（setosa, versicolor, virginica）之间平均分配，这些数据在数据集中是相邻的。你的训练数据包含了所有的setosa和versicolor，而测试集则包含了所有的virginica数据。没有重叠，这会导致样本外的问题。最好平衡你的训练和测试集以避免这种情况。- 学习率和模型数量的组合对我来说看起来太低了。拟合收敛为
(1-lr)^n。当lr = 1e-3和n = 1000时，你只能模拟63.2%的数据幅度。也就是说，即使每个模型对每个样本的预测都正确，你也只能估计出63.2%的正确值。用平均值而不是0来初始化拟合会有所帮助，因为这样效果是回归到均值，而不是仅仅拖累拟合值。
