R: implémenter mon propre algorithme de renforcement de gradient

J'essaie d'écrire mon propre algorithme de renforcement de gradient. Je comprends qu'il existe des packages existants comme gbmet xgboost,mais je voulais comprendre comment l'algorithme fonctionne en écrivant le mien.

J'utilise l' irisensemble de données et mon résultat est Sepal.Length(continu). Ma fonction de perte est mean(1/2*(y-yhat)^2)(essentiellement l'erreur quadratique moyenne avec 1/2 devant), donc mon gradient correspondant n'est que le résiduel y - yhat. J'initialise les prédictions à 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)
}

Avec cela, j'ai divisé l' irisensemble de données en un ensemble de données de formation et de test et j'adapte mon modèle à celui-ci.

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)

Maintenant, je calcule les valeurs prévues à partir de my.model. Pour my.model, les valeurs ajustées sont 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

J'ai quelques questions

1. Mon algorithme d'amplification du gradient semble-t-il correct?
2. Ai-je calculé yhats.mymodcorrectement les valeurs prévues ?

1. Oui, cela semble correct. À chaque étape, vous ajustez les pseudo-résidus, qui sont calculés comme la dérivée de la perte par rapport à l'ajustement. Vous avez correctement dérivé ce gradient au début de votre question, et même pris la peine d'obtenir le bon facteur 2.
2. Cela semble également correct. Vous agrégez les modèles, pondérés par le taux d'apprentissage, comme vous l'avez fait pendant la formation.

Mais pour répondre à quelque chose qui n'a pas été demandé, j'ai remarqué que votre configuration d'entraînement présente quelques bizarreries.

• L' irisensemble de données est réparti également entre 3 espèces (setosa, versicolor, virginica) et celles-ci sont adjacentes dans les données. Vos données d'entraînement ont tous les setosa et versicolor, tandis que l'ensemble de test a tous les exemples virginica. Il n'y a pas de chevauchement, ce qui entraînera des problèmes hors échantillon. Il est préférable d'équilibrer votre entraînement et vos tests pour éviter cela.
• La combinaison du taux d'apprentissage et du nombre de modèles me semble trop faible. L'ajustement converge au fur et à mesure (1-lr)^n. Avec lr = 1e-3et n = 1000vous ne pouvez modéliser que 63,2% de la magnitude des données. Autrement dit, même si chaque modèle prédit correctement chaque échantillon, vous estimeriez 63,2% de la valeur correcte. L'initialisation de l'ajustement avec une moyenne, au lieu de 0, serait utile car alors l'effet est une régression vers la moyenne au lieu d'une simple traînée.