R: implementar mi propio algoritmo de aumento de gradiente

Estoy tratando de escribir mi propio algoritmo de aumento de gradiente. Entiendo que hay paquetes existentes como gbmy, xgboost,pero quería entender cómo funciona el algoritmo escribiendo el mío.

Estoy usando el irisconjunto de datos y mi resultado es Sepal.Length(continuo). Mi función de pérdida es mean(1/2*(y-yhat)^2)(básicamente el error cuadrático medio con 1/2 en el frente), por lo que mi gradiente correspondiente es solo el residual y - yhat. Estoy inicializando las predicciones en 0.

library(rpart)
data(iris)

#Define gradient
grad.fun <- function(y, yhat) {return(y - yhat)}

mod <- list()

grad_boost <- function(data, learning.rate, M, grad.fun) {
  # Initialize fit to be 0
  fit <- rep(0, nrow(data))
  grad <- grad.fun(y = data$Sepal.Length, yhat = fit)

  # Initialize model
  mod[[1]] <- fit

  # Loop over a total of M iterations
  for(i in 1:M){

    # Fit base learner (tree) to the gradient
    tmp <- data$Sepal.Length
    data$Sepal.Length <- grad
    base_learner <- rpart(Sepal.Length ~ ., data = data, control = ("maxdepth = 2"))
    data$Sepal.Length <- tmp

    # Fitted values by fitting current model
    fit <- fit + learning.rate * as.vector(predict(base_learner, newdata = data))

    # Update gradient
    grad <- grad.fun(y = data$Sepal.Length, yhat = fit)

    # Store current model (index is i + 1 because i = 1 contain the initialized estiamtes)
    mod[[i + 1]] <- base_learner

  }
  return(mod)
}

Con esto, dividí el irisconjunto de datos en un conjunto de datos de entrenamiento y prueba y ajusté mi modelo a él.

train.dat <- iris[1:100, ]
test.dat <- iris[101:150, ]
learning.rate <- 0.001
M = 1000
my.model <- grad_boost(data = train.dat, learning.rate = learning.rate, M = M, grad.fun = grad.fun)

Ahora calculo los valores predichos de my.model. Para my.model, los valores ajustados son 0 (vector of initial estimates) + learning.rate * predictions from tree 1 + learning rate * predictions from tree 2 + ... + learning.rate * predictions from tree M.

yhats.mymod <- apply(sapply(2:length(my.model), function(x) learning.rate * predict(my.model[[x]], newdata = test.dat)), 1, sum)

# Calculate RMSE
> sqrt(mean((test.dat$Sepal.Length - yhats.mymod)^2))
[1] 2.612972

Tengo algunas preguntas

1. ¿Mi algoritmo de aumento de gradiente se ve bien?
2. ¿Calculé los valores pronosticados yhats.mymodcorrectamente?

1. Sí, esto parece correcto. En cada paso, se ajusta a los residuos psuedo, que se calculan como la derivada de la pérdida con respecto al ajuste. Has derivado correctamente este gradiente al comienzo de tu pregunta, e incluso te has molestado en obtener el factor 2 correcto.
2. Esto también se ve correcto. Se está agregando a través de los modelos, ponderado por la tasa de aprendizaje, tal como lo hizo durante el entrenamiento.

Pero para abordar algo que no se le preguntó, noté que su configuración de entrenamiento tiene algunas peculiaridades.

• El irisconjunto de datos se divide por igual entre 3 especies (setosa, versicolor, virginica) y estos son adyacentes en los datos. Sus datos de entrenamiento tienen todos los setosa y versicolor, mientras que el conjunto de prueba tiene todos los ejemplos de virginica. No hay superposición, lo que dará lugar a problemas fuera de la muestra. Es preferible equilibrar su entrenamiento y pruebas para evitar esto.
• La combinación de tasa de aprendizaje y conteo de modelos me parece demasiado baja. El ajuste converge como (1-lr)^n. Con lr = 1e-3y n = 1000solo puede modelar el 63.2% de la magnitud de los datos. Es decir, incluso si cada modelo predice cada muestra correctamente, estaría estimando el 63.2% del valor correcto. Inicializar el ajuste con un promedio, en lugar de 0, ayudaría desde entonces, el efecto es una regresión a la media en lugar de solo un arrastre.