Algoritmo genético para selección de predictores
Joaquín Amat Rodrigo j.amatrodrigo@gmail.com
Noviembre, 2018 actualizado Enero, 2020
Versión PDF: Github
Más sobre ciencia de datos: joaquinamatrodrigo.github.io o cienciadedatos.net
Introducción
Los algoritmos genéticos son métodos de optimización heurística que, entre otras aplicaciones, pueden emplearse para encontrar la combinación de variables que consigue maximizar la capacidad predictiva de un modelo. Su funcionamiento está inspirado en la teoría evolutiva de selección natural propuesta por Darwin y Alfred Russel: los individuos de una población se reproducen generando nuevos descendientes, cuyas características, son combinación de las características de los progenitores (más ciertas mutaciones). De todos ellos, únicamente los mejores individuos sobreviven y pueden reproducirse de nuevo, transmitiendo así sus características a las siguientes generaciones.
Los algoritmos genéticos son solo una de las muchas estrategias que existen para seleccionar los predictores más relevantes, y no tiene por qué ser la más adecuada en todos los escenarios. Por ejemplo, existen estrategias iterativas Stepwise selection, modelos como Random Forest, Boosting y Lasso capaces de excluir automáticamente predictores, y técnicas de reducción de dimensión como PCA y t-SNE.
La implementación de algoritmo genético que se muestra en este documento pretende ser lo más explicativa posible aunque para ello no sea la más eficiente.
El código de las funciones desarrolladas a lo largo del documento puede descargarse en el siguiente link.
Algoritmo
Aunque existen variaciones, algunas de las cuales se describen a lo largo de este documento, en términos generales, la estructura de un algoritmo genético para la selección de predictores sigue los siguientes pasos:
Crear una población inicial aleatoria de \(P\) individuos. En este caso, cada individuo representa una combinación de predictores.
Calcular la fortaleza (fitness) de cada individuo de la población.
Crear una nueva población vacía y repetir los siguientes pasos hasta que se hayan creado \(P\) nuevos individuos.
3.1 Seleccionar dos individuos de la población existente, donde la probabilidad de selección es proporcional al fitness de los individuos.
3.2 Cruzar los dos individuos seleccionados para generar un nuevo descendiente (crossover).
3.3 Aplicar un proceso de mutación aleatorio sobre el nuevo individuo.
3.4 Añadir el nuevo individuo a la nueva población.
Reemplazar la antigua población por la nueva.
Si no se cumple un criterio de parada, volver al paso 2.
En los siguientes apartados se describe cada una de las etapas del proceso para, finalmente, combinarlas todas en una única función.
Población inicial
En el contexto de algoritmos genéticos, el término individuo hace referencia a cada una de las posibles soluciones del problema que se quiere optimizar. En el caso particular de la selección de mejores predictores, cada individuo representa una posible combinación de variables. Para representar dichas combinaciones, se pueden emplear vectores binarios, cuya longitud es igual al número total de predictores disponibles, y cada posición toma el valor TRUE/FALSE dependiendo de si el predictor que ocupa esa posición se incluye o excluye. También es común encontrar la representación en términos 0/1.
Por ejemplo, supóngase que las variables disponibles son: \(X_1\), \(X_2\), \(X_3\), \(X_4\) y \(X_5\). El individuo \(TRUE, FALSE, TRUE, TRUE, FALSE\), o su equivalente codificación \(1, 0, 1, 1, 0\) representa la selección de los predictores \(X_1\), \(X_3\), y \(X_4\).
El primer paso del algoritmo genético para la selección de predictores consiste en crear una población inicial aleatoria de individuos. La siguiente función crea una matriz binaria en la que, cada fila, está formada una combinación aleatoria de valores TRUE y FALSE. Además, el número máximo y mínimo de TRUEs por fila puede estar acotado. Esta acotación resulta útil cuando se quiere limitar, en cierta medida, el número de predictores que pueden incluir los individuos. Es en cierta medida porque, debido a los cruces y mutaciones a lo largo de las generaciones, se pueden crear individuos que excedan los límites iniciales.
crear_poblacion <- function(n_poblacion,
n_variables,
n_max = NULL,
n_min = NULL,
verbose = TRUE) {
# Parameters
# ----------
# n_poblacion: "integer"
# número total de individuos de la población.
#
# n_variables: "integer"
# longitud de los individuos.
#
# n_max: "integer"
# número máximo de TRUEs que puede contener un individuo.
#
# n_min: "integer"
# número mínimo de TRUEs que puede contener un individuo.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# "matrix" de tamaño n_poblacion x n_variables que representa una población.
# COMPROBACIONES
# ----------------------------------------------------------------------------
if (isTRUE(n_max > n_variables)) {
stop("n_max no puede ser mayor que n_variables.")
}
# Si no se especifica n_max, el número máximo de predictores (TRUEs) que puede
# contener un individuo es igual al número total de variables disponibles.
if (is.null(n_max)) {
n_max <- n_variables
}
# Si no se especifica n_min, el número mínimo de predictores (TRUEs) que puede
# contener un individuo es 1.
if (is.null(n_min)) {
n_min <- 1
}
# CREAR POBLACIÓN
# ----------------------------------------------------------------------------
# Matriz donde almacenar los individuos generados.
poblacion <- matrix(data = NA, nrow = n_poblacion, ncol = n_variables)
# Bucle para crear cada individuo.
for (i in 1:n_poblacion) {
# Se selecciona (con igual probabilidad) el número de valores TRUE que puede
# tener el individuo, dentro del rango acotado por n_min y n_max.
n_true <- sample(x = n_min:n_max, size = 1)
# Se crea un vector con todo FALSE que representa el individuo.
individuo <- rep(FALSE, times = n_variables)
# Se sustituyen n_true posiciones aleatorias por valores TRUE.
individuo[sample(x = 1:n_variables, size = n_true)] <- TRUE
# Se añade el nuevo individuo a la población.
poblacion[i, ] <- individuo
}
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(poblacion, 'fecha_creacion') <- Sys.time()
attr(poblacion, 'numero_individuos') <- n_poblacion
attr(poblacion, "class") <- c(attr(poblacion, "class"), "poblacion")
if (verbose) {
cat("Población inicial creada", "\n")
cat("------------------------", "\n")
cat("Fecha creación:", as.character(Sys.time()), "\n")
cat("Número de individuos =", n_poblacion, "\n")
cat("Número de predictores mínimo por individuo =", n_min, "\n")
cat("Número de predictores máximo por individuo =", n_max, "\n")
cat("\n")
}
return(poblacion)
}Ejemplo
Se crea una población de 10 individuos de longitud 8, con un número de valores TRUE acotado entre 1 y 5.
poblacion <- crear_poblacion(
n_poblacion = 10,
n_variables = 8,
n_max = 5,
n_min = 1,
verbose = TRUE
)## Población inicial creada
## ------------------------
## Fecha creación: 2020-02-01 18:58:09
## Número de individuos = 10
## Número de predictores mínimo por individuo = 1
## Número de predictores máximo por individuo = 5Fitness de un individuo
Cada individuo de la población debe ser evaluado para cuantificar su fortaleza (fitness). Dado que, en este caso, el objetivo es encontrar la combinación de predictores que da lugar al mejor modelo, el fitness de un individuo se calcula con una métrica de calidad del modelo. Dependiendo de la métrica, la relación con el fitness puede ser:
Error del modelo: si se emplea una métrica de error, el individuo tiene mayor fitness cuanto menor sea el error.
Accuracy, precision, recall, F1…: con este tipo de métricas, el individuo tiene mayor fitness cuanto mayor sea la métrica.
Para conseguir que, independientemente de la métrica, cuanto mayor sea su valor, mayor el fitness del individuo, se puede utilizar el \(-(error)\). De esta forma, la estrategia de selección es la misma. Es importante recordar que, en el caso del \(-(error)\), el intervalo de valores posibles es \(-\infty\) a 0. Además, para que las estimaciones sean robustas y evitar el overfitting, es muy importante recurrir a estrategias de validación cruzada o bootstrapping en la cuantificación del fitness.
En este ejemplo, se presentan tres opciones de modelos de evaluación: un modelo lineal por mínimos cuadrados, un modelo de regresión logística y un random forest.
El modelo lineal solo puede utilizarse para problemas de regresión y emplea la métrica -MSE como indicativo de fitness.
El modelo de regresión logística solo puede utilizarse para clasificación binaria y emplea las métricas Accuracy, índice Kappa o índice F1.
El random forest puede utilizarse para problemas de regresión, en cuyo caso emplea la métrica -(MSE), y clasificación, en cuyo caso emplea el Accuracy, índice Kappa o índice F1.
Versión con random forest
calcular_fitness_individuo_rf <- function(x,
y,
cv,
metrica = NULL,
nivel_referencia = NULL,
n_tree = 50,
seed = 123,
verbose = TRUE,
...) {
# Parameters
# ----------
# x: "matrix", "data.frame", tibble
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# metrica: ("-mse", -mae, "f1", "kappa", "accuracy")
# métrica utilizada para calcular el fitness. Por defecto es el -mse para
# regresión y accuracy para clasificación. En el caso de clasificación
# binaria, se acepta también f1 y kappa.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# n_tree: "integer"
# número de árboles del modelo randomforest.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# fitness ("numeric") del individuo obtenido por validación cruzada empleando
# random forest como modelo.
# LIBRERÍAS
# ----------------------------------------------------------------------------
library(MLmetrics)
library(caret)
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (!is.null(metrica)) {
if (metrica %in% c("-mse", "-mae") && !is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de regresión.")
)
}
if (metrica %in% c("accuracy", "f1", "kappa") && is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de clasificación.")
)
}
if (metrica %in% c("f1", "kappa") && length(unique(y)) > 2) {
stop(
paste(
"La métrica",
metrica,
"no es válida para problemas de clasificación con más de dos clases."
)
)
}
if (metrica %in% c("f1", "kappa") && is.null(nivel_referencia)) {
stop(
"Para la métrica kappa es necesario indicar el nivel de referencia."
)
}
}
if (is.null(metrica)) {
if (is.numeric(y)) {
metrica <- "-mse"
}else{
metrica <- "accuracy"
}
}
# REPARTO DE LAS OBSERVACIONES PARA LA VALIDACIÓN CRUZADA (CV)
# ----------------------------------------------------------------------------
set.seed(seed)
indices_cv <- caret::createFolds(y = y, k = cv, list = FALSE)
# BUCLE DE ENTRENAMIENTO Y VALIDACIÓN CRUZADA (CV)
# ----------------------------------------------------------------------------
# Vector para almacenar el fitness de cada iteración CV.
metrica_cv <- rep(NA, times = cv)
for (i in 1:cv) {
set.seed(seed)
# modelo <- randomForest::randomForest(
# x = x[indices_cv != i, , drop = FALSE],
# y = y[indices_cv != i],
# ntree = n_tree,
# importance = FALSE,
# localImp = FALSE
# )
# predicciones <- predict(modelo, data = x[indices_cv == i, , drop = FALSE])
modelo <- ranger::ranger(
x = as.data.frame(x[indices_cv != i, , drop = FALSE]),
y = y[indices_cv != i],
num.trees = n_tree,
importance = "none",
keep.inbag = FALSE,
oob.error = FALSE
)
predicciones <- predict(
modelo,
data = as.data.frame(x[indices_cv == i, , drop = FALSE])
)
predicciones <- predicciones$predictions
if (metrica == "-mse") {
residuos <- predicciones - y[indices_cv == i]
metrica_cv[i] <- -(mean(residuos^2))
} else if (metrica == "-mae") {
residuos <- predicciones - y[indices_cv == i]
metrica_cv[i] <- -(mean(abs(residuos)))
} else if (metrica == "accuracy") {
accuracy <- MLmetrics::Accuracy(
y_pred = predicciones,
y_true = y[indices_cv == i]
)
metrica_cv[i] <- accuracy
} else if (metrica == "kappa") {
mat_confusion <- caret::confusionMatrix(
table(
predicciones,
y[indices_cv == i]
),
positive = nivel_referencia
)
kappa <- mat_confusion$overall["Kappa"]
metrica_cv[i] <- kappa
} else if (metrica == "f1") {
if (sum(predicciones == nivel_referencia) < 1) {
print(
paste(
"Ningún valor predicho ha sido el nivel de referencia,",
"por lo tanto, el valor de precisión no puede ser calculado",
"y tampoco la métrica F1."
)
)
} else {
f1 <- MLmetrics::F1_Score(
y_true = y[indices_cv == i],
y_pred = predicciones,
positive = nivel_referencia
)
metrica_cv[i] <- f1
}
}
# Información detalla de la partición en caso de que el cálculo de la
# métrica falle.
if (verbose) {
if (is.na(metrica_cv[i])) {
print(paste("En la particion", i, "de CV, fitness del individuo no",
"ha podido ser calculado."
))
print("Valores reales:")
print(y[indices_cv == i])
print("Valores predichos:")
print(predicciones)
}
}
}
# COMPROBACIÓN DE LAS MÉTRICAS OBTENIDAS POR CV
# ------------------------------------------------------------------------
if (any(is.na(metrica_cv))) {
warning(paste(
"En una o más particiones de CV, la métrica del individuo no",
"ha podido ser calculado."
))
}
if (is.na(mean(metrica_cv, na.rm = TRUE))) {
stop("El fitness del individuo no ha podido ser calculado.")
}
valor_metrica <- mean(metrica_cv, na.rm = TRUE)
fitness <- mean(metrica_cv, na.rm = TRUE)
modelo <- "Random Forest"
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ------------------------------------------------------------------------
if (verbose) {
cat("El individuo ha sido evaluado", "\n")
cat("-----------------------------", "\n")
cat("Métrica:", metrica, "\n")
cat("Valor métrica:", valor_metrica, "\n")
cat("Fitness:", fitness, "\n")
cat("Modelo de evaluación:", modelo, "\n")
cat("\n")
}
return(fitness)
}Versión con modelo lineal
calcular_fitness_individuo_lm <- function(x,
y,
cv,
metrica = NULL,
seed = 123,
verbose = TRUE,
...) {
# Parameters
# ----------
# x: "matrix", "data.frame", "tibble"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# metrica: ("-mse", "-mae)
# métrica utilizada para calcular el fitness. Por defecto se emplea -mse.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# fitness ("numeric") del individuo obtenido por validación cruzada empleando
# regresión lineal por mínimos cuadrados modelo.
# LIBRERÍAS
# ------------------------------------------------------------------------
library(MLmetrics)
library(caret)
# COMPROBACIONES INICIALES
# ------------------------------------------------------------------------
if (!is.numeric(y)) {
stop(
"El método lm solo es valido para problemas de regresión."
)
}
if (!is.null(metrica)) {
if (!metrica %in% c("-mse", "-mae")) {
stop(
paste(
"La métrica", metrica, "no es válida para problemas de regresión.",
"Emplear el -mse o -mae."
)
)
}
}
if (is.null(metrica)) {
metrica <- "-mse"
}
# REPARTO DE LAS OBSERVACIONES PARA LA VALIDACIÓN CRUZADA (CV)
# ------------------------------------------------------------------------
set.seed(seed)
indices_cv <- caret::createFolds(y = y, k = cv, list = FALSE)
# BUCLE DE ENTRENAMIENTO Y VALIDACIÓN CRUZADA (CV)
# ------------------------------------------------------------------------
# Vector para almacenar la métrica de cada iteración CV.
metrica_cv <- rep(NA, times = cv)
for (i in 1:cv) {
datos_modelo <- as.data.frame(
cbind(
x[indices_cv != i, , drop = FALSE],
y = y[indices_cv != i]
)
)
modelo <- lm(
formula = y ~ .,
data = datos_modelo
)
predicciones <- predict(
modelo,
newdata = as.data.frame(
x[indices_cv == i, , drop = FALSE]
)
)
residuos <- predicciones - y[indices_cv == i]
if (metrica == "-mse") {
metrica_cv[i] <- -(mean(residuos^2))
}
if (metrica == "-mae") {
metrica_cv[i] <- -(mean(abs(residuos)))
}
# Información detalla de la partición en caso de que el cálculo de la
# métrica falle.
if (verbose) {
if (is.na(metrica_cv[i])) {
print(paste("En la particion", i, "de CV, fitness del individuo no",
"ha podido ser calculado."
))
print("Valores reales:")
print(y[indices_cv == i])
print("Valores predichos:")
print(predicciones)
}
}
}
# COMPROBACIÓN DE LAS MÉTRICAS OBTENIDAS POR CV
# ------------------------------------------------------------------------
if (any(is.na(metrica_cv))) {
warning(paste(
"En una o más particiones de CV, la métrica del individuo no",
"ha podido ser calculado."
))
}
if (is.na(mean(metrica_cv, na.rm = TRUE))) {
stop("El fitness del individuo no ha podido ser calculado.")
}
valor_metrica <- mean(metrica_cv, na.rm = TRUE)
fitness <- mean(metrica_cv, na.rm = TRUE)
modelo <- "Regresión lineal (lm)"
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ------------------------------------------------------------------------
if (verbose) {
cat("El individuo ha sido evaluado", "\n")
cat("-----------------------------", "\n")
cat("Métrica:", metrica, "\n")
cat("Valor métrica:", valor_metrica, "\n")
cat("Fitness:", fitness, "\n")
cat("Modelo de evaluación:", modelo, "\n")
cat("\n")
}
return(fitness)
}Versión con modelo de regresión logística
calcular_fitness_individuo_glm <- function(x,
y,
cv,
metrica = NULL,
nivel_referencia = NULL,
seed = 123,
verbose = TRUE,
...) {
# Parameters
# ----------
# x: "matrix", "data.frame" , "tibble"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# metrica: {f1, kappa, accuracy}
# métrica utilizada para calcular el fitness. Por defecto es accuracy.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# fitness ("numeric") del individuo obtenido por validación cruzada empleando
# regresión logística como modelo.
# LIBRERÍAS
# ------------------------------------------------------------------------
library(MLmetrics)
library(caret)
# COMPROBACIONES INICIALES
# ------------------------------------------------------------------------
if (!is.null(metrica)) {
if (!metrica %in% c("accuracy", "f1", "kappa")) {
stop(
"Las métricas permitidas con el modelo glm son: accuracy, f1 y kappa."
)
}
}
if (is.null(metrica)) {
metrica <- "accuracy"
}
if (!is.character(y) & !is.factor(y)) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación,",
"la variable respuesta debe ser character o factor."
)
)
}
if (is.character(y)) {
y <- as.factor(y)
}
if (length(levels(y)) != 2) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación binaria,",
"la variable respuesta debe tener dos niveles."
)
)
}
# REPARTO DE LAS OBSERVACIONES PARA LA VALIDACIÓN CRUZADA (CV)
# ------------------------------------------------------------------------
set.seed(seed)
indices_cv <- caret::createFolds(y = y, k = cv, list = FALSE)
# BUCLE DE ENTRENAMIENTO Y VALIDACIÓN CRUZADA (CV)
# ------------------------------------------------------------------------
# Vector para almacenar la métrica de cada iteración CV.
metrica_cv <- rep(NA, times = cv)
for (i in 1:cv) {
datos_modelo <- as.data.frame(
cbind(
x[indices_cv != i, , drop = FALSE],
y = y[indices_cv != i]
)
)
modelo <- glm(
formula = y ~ .,
family = "binomial",
data = datos_modelo
)
predicciones <- predict(
modelo,
newdata = as.data.frame(
x[indices_cv == i, , drop = FALSE]
),
type = "response"
)
predicciones <- unname(predicciones)
predicciones <- ifelse(predicciones < 0.5, levels(y)[1], levels(y)[2])
predicciones <- factor(x = predicciones, levels = levels(y))
if (metrica == "accuracy") {
accuracy <- MLmetrics::Accuracy(
y_pred = predicciones,
y_true = y[indices_cv == i]
)
metrica_cv[i] <- accuracy
} else if (metrica == "kappa") {
mat_confusion <- caret::confusionMatrix(
table(
predicciones,
y[indices_cv == i]
),
positive = nivel_referencia
)
kappa <- mat_confusion$overall["Kappa"]
metrica_cv[i] <- kappa
} else if (metrica == "f1") {
if (sum(predicciones == nivel_referencia) < 1) {
print(
paste(
"Ningún valor predicho ha sido el nivel de referencia,",
"por lo tanto, el valor de precisión no puede ser calculado",
"y tampoco la métrica F1."
)
)
} else {
f1 <- MLmetrics::F1_Score(
y_true = y[indices_cv == i],
y_pred = predicciones,
positive = nivel_referencia
)
metrica_cv[i] <- f1
}
}
# Información detalla de la partición en caso de que el cálculo de la
# métrica falle.
if (verbose) {
if (is.na(metrica_cv[i])) {
print(paste("En la particion", i, "de CV, fitness del individuo no",
"ha podido ser calculado."
))
print("Valores reales:")
print(y[indices_cv == i])
print("Valores predichos:")
print(predicciones)
}
}
}
# COMPROBACIÓN DE LAS MÉTRICAS OBTENIDAS POR CV
# ------------------------------------------------------------------------
if (any(is.na(metrica_cv))) {
warning(paste(
"En una o más particiones de CV, la métrica del individuo no",
"ha podido ser calculado."
))
}
if (is.na(mean(metrica_cv, na.rm = TRUE))) {
stop("El fitness del individuo no ha podido ser calculado.")
}
valor_metrica <- mean(metrica_cv, na.rm = TRUE)
fitness <- mean(metrica_cv, na.rm = TRUE)
modelo <- "Regresión Logística (glm)"
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ------------------------------------------------------------------------
if (verbose) {
cat("El individuo ha sido evaluado", "\n")
cat("-----------------------------", "\n")
cat("Métrica:", metrica, "\n")
cat("Valor métrica:", valor_metrica, "\n")
cat("Fitness:", fitness, "\n")
cat("Modelo de evaluación:", modelo, "\n")
cat("\n")
}
return(fitness)
}Ejemplo
Se calcula el fitness de un individuo que incluye 10 predictores simulados.
set.seed(123)
predictores <- matrix(data = rnorm(n = 1000), nrow = 100, ncol = 10)
var_respuesta <- rnorm(n = 100)
calcular_fitness_individuo_rf(
x = predictores,
y = var_respuesta,
metrica = "-mse",
cv = 5,
seed = 123,
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -1.211097
## Fitness: -1.211097
## Modelo de evaluación: Random Forest## [1] -1.211097calcular_fitness_individuo_lm(
x = predictores,
y = var_respuesta,
metrica = "-mse",
cv = 5,
seed = 123,
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -1.244765
## Fitness: -1.244765
## Modelo de evaluación: Regresión lineal (lm)## [1] -1.244765Fitness de todos los individuos de una población
Esta función recibe como argumentos una población de individuos, una matriz de predictores, un vector con la variable respuesta, y devuelve el fitness de cada individuo.
calcular_fitness_poblacion <- function(poblacion,
x,
y,
modelo,
cv,
seed = 123,
metrica = NULL,
nivel_referencia = NULL,
verbose = TRUE) {
# Parameters
# ----------
# poblacion: "matrix"
# matriz boleana que representa la población de individuos.
#
# x: "matrix" o "data.frame"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# modelo: ("lm", "glm", "randomforest").
# tipo de modelo empleado para calcular el fitness.
#
# metrica: ("-mse", "-mae", "f1", "kappa", "accuracy").
# métrica utilizada para calcular el fitness. Por defecto es el -mse para
# regresión y accuracy para clasificación. En el caso de clasificación
# binaria, se acepta también f1 y kappa.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
#
# Return
# ----------
# Vector ("numeric") con el fitness de todos los individuos de la población,
# obtenido por validación cruzada. El orden de los valores se corresponde con
# el orden de las filas de la matriz población.
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (!modelo %in% c("lm", "glm", "randomforest")) {
stop(
paste("El modelo empleado para calcular el fitness debe ser",
"lm (regresion lineal)",
"glm (regresión logistica)",
"o randomforest."
)
)
}
if (ncol(poblacion) != ncol(x)) {
stop(
"n_variables de la población debe ser igual al número de columnas de x."
)
}
# Si el modelo es de regresión lineal
if (modelo == "lm") {
if (!is.numeric(y)) {
stop(
"El método lm solo es valido para problemas de regresión."
)
}
if (!is.null(metrica)) {
if (!metrica %in% c("-mse", "-mae")) {
stop(
paste(
"La métrica", metrica, "no es válida para problemas de regresión.",
"Emplear el -mse o -mae."
)
)
}
}else{
# Por defecto se emplea -mse
metrica <- "-mse"
}
}
# Si el modelo es de regresión logística
if (modelo == "glm") {
if (!is.null(metrica)) {
if (!metrica %in% c("accuracy", "f1", "kappa")) {
stop(
"Las métricas permitidas con el modelo glm son: accuracy, f1 y kappa."
)
}
}else{
# Por defecto se emplea accuracy
metrica <- "accuracy"
}
if (!is.character(y) & !is.factor(y)) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación,",
"la variable respuesta debe ser character o factor."
)
)
}
if (is.character(y)) {
y <- as.factor(y)
}
if (length(levels(y)) != 2) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación binaria,",
"la variable respuesta debe tener dos niveles."
)
)
}
}
# Si el modelo es de random forest
if (modelo == "randomforest") {
if (!is.null(metrica)) {
if (metrica %in% c("-mse", "-mae") && !is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de regresión.")
)
}
if (metrica %in% c("accuracy", "f1", "kappa") && is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de clasificación.")
)
}
if (metrica %in% c("f1", "kappa") && length(unique(y)) > 2) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de clasificación,",
"binaria."
)
)
}
if (metrica %in% c("f1", "kappa") && is.null(nivel_referencia)) {
stop(
"Para la métrica kappa es necesario indicar el nivel de referencia."
)
}
}
if (is.null(metrica)) {
if (is.numeric(y)) {
metrica <- "-mse"
}else{
metrica <- "accuracy"
}
}
}
# SELECCIÓN FUNCIÓN AUXILIAR PARA CALCULAR EL FITNESS
# ----------------------------------------------------------------------------
# Tipo de modelo utilizado para calcular el fitness.
if (modelo == "lm") {
calcular_fitness_individuo <- calcular_fitness_individuo_lm
} else if (modelo == "randomforest") {
calcular_fitness_individuo <- calcular_fitness_individuo_rf
} else if (modelo == "glm") {
calcular_fitness_individuo <- calcular_fitness_individuo_glm
}
# CÁLCULO DEL FITNNES DE CADA INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
# Vector donde almacenar el fitness de cada individuo.
fitness_poblacion <- rep(NA, times = nrow(poblacion))
for (i in 1:nrow(poblacion)) {
individuo <- poblacion[i, ]
if (verbose) {
cat("Individuo", i, ":", individuo, "\n")
}
fitness_individuo <- calcular_fitness_individuo(
x = x[, individuo, drop = FALSE],
y = y,
cv = cv,
metrica = metrica,
nivel_referencia = nivel_referencia,
seed = seed,
verbose = verbose
)
fitness_poblacion[i] <- fitness_individuo
}
# MEJOR INDIVIDUO DE LA POBLACIÓN
# ------------------------------------------------------------------------
# Se identifica el mejor individuo de toda la población, el de mayor
# fitness.
indice_mejor_individuo <- which.max(fitness_poblacion)
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("------------------", "\n")
cat("Población evaluada", "\n")
cat("------------------", "\n")
cat("Modelo empleado para el cálculo del fitness:", modelo, "\n")
cat("Métrica empleado para el cálculo del fitness:", metrica, "\n")
cat("Mejor fitness encontrado:", fitness_poblacion[indice_mejor_individuo], "\n")
cat("Mejor secuencia encontrada:", "\n")
print(poblacion[indice_mejor_individuo,])
cat("Mejores predictores encontrados:", "\n")
print(colnames(x)[poblacion[indice_mejor_individuo,]])
cat("\n")
}
return(fitness_poblacion)
}Ejemplo
Se calcula el fitness de todos los individuos de una población formada por 10 individuos de longitud 8, con un número de valores TRUE acotado entre 1 y 5.
# Población simulada
poblacion <- crear_poblacion(
n_poblacion = 10,
n_variables = 8,
n_max = 5,
n_min = 1,
verbose = FALSE
)
# Datos simulados
set.seed = 123
predictores <- as.data.frame(matrix(data = rnorm(n = 1000), nrow = 100, ncol = 8))
var_respuesta <- rnorm(n = 100)
# Cálculo del fitness de todos los individuos
fitness_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
x = predictores,
y = var_respuesta,
modelo = "lm",
cv = 5,
seed = 123,
verbose = TRUE
)## Individuo 1 : FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9723512
## Fitness: -0.9723512
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 2 : FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9719444
## Fitness: -0.9719444
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 3 : TRUE FALSE TRUE TRUE FALSE TRUE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9850946
## Fitness: -0.9850946
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 4 : TRUE FALSE FALSE FALSE FALSE TRUE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9923327
## Fitness: -0.9923327
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 5 : FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9979785
## Fitness: -0.9979785
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 6 : FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9720212
## Fitness: -0.9720212
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 7 : FALSE TRUE TRUE TRUE FALSE TRUE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9834769
## Fitness: -0.9834769
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 8 : TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -1.022854
## Fitness: -1.022854
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 9 : FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9822012
## Fitness: -0.9822012
## Modelo de evaluación: Regresión lineal (lm)
##
## Individuo 10 : FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -0.9902582
## Fitness: -0.9902582
## Modelo de evaluación: Regresión lineal (lm)
##
## ------------------
## Población evaluada
## ------------------
## Modelo empleado para el cálculo del fitness: lm
## Métrica empleado para el cálculo del fitness: -mse
## Mejor fitness encontrado: -0.9719444
## Mejor secuencia encontrada:
## [1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE
## Mejores predictores encontrados:
## [1] "V6" "V7"fitness_poblacion## [1] -0.9723512 -0.9719444 -0.9850946 -0.9923327 -0.9979785 -0.9720212
## [7] -0.9834769 -1.0228537 -0.9822012 -0.9902582El vector devuelto contiene el fitness de cada uno de los individuos en el mismo orden que se encuentran en la matriz de la población.
Seleccionar individuos
La forma en que se seleccionan los individuos que participan en cada cruce difiere en las distintas implementaciones de los algoritmos genéticos. Por lo general, todas ellas tienden a favorecer la selección de aquellos individuos con mayor fitness. Algunas de las estrategias más comunes son:
Método de ruleta: la probabilidad de que un individuo sea seleccionado es proporcional a su fitness relativo, es decir, a su fitness dividido por la suma del fitness de todos los individuos de la población. Si el fitness de un individuo es el doble que el de otro, también lo será la probabilidad de que sea seleccionado. Este método presenta problemas si el fitness de unos pocos individuos es muy superior (varios órdenes de magnitud) al resto, ya que estos serán seleccionados de forma repetida y casi todos los individuos de la siguiente generación serán “hijos” de los mismos “padres” (poca variación).
Método rank: la probabilidad de selección de un individuo es inversamente proporcional a la posición que ocupa tras ordenar todos los individuos de mayor a menor fitness. Este método es menos agresivo que el método ruleta cuando la diferencia entre los mayores fitness es varios órdenes de magnitud superior al resto.
Selección competitiva (tournament): se seleccionan aleatoriamente dos parejas de individuos de la población (todos con la misma probabilidad). De cada pareja se selecciona el que tenga mayor fitness. Finalmente, se comparan los dos finalistas y se selecciona el de mayor fitness. Este método tiende a generar una distribución de la probabilidad de selección más equilibrada que las dos anteriores.
Selección truncada (truncated selection): se realizan selecciones aleatorias de individuos, habiendo descartado primero los n individuos con menor fitness de la población.
La siguiente función recibe como argumento un vector con el fitness de cada individuo y selecciona una de las posiciones, donde la probabilidad de selección es proporcional al fitness relativo.
La conversión de fitness a probabilidad es distinta dependiendo de la métrica utilizada para calcular el fitness.
Si el fitness es el valor negativo del error (-MSE), cuanto menor sea el fitness (menor la magnitud del valor negativo), mayor debe ser la probabilidad de ser seleccionado. Para lograr la conversión se emplea \(fitness = \frac{-1}{-1*(MSE)}\).
Si el fitness es el accuracy, f1 o kappa, cuanto mayor sea el fitness, mayor debe ser la probabilidad de ser seleccionado.
seleccionar_individuo <- function(vector_fitness,
metrica,
metodo_seleccion = "tournament",
verbose = TRUE) {
# Parameters
# ----------
# vector_fitness: "numeric"
# un vector con el fitness de cada individuo.
#
# metrica: ("-mse", "-mae", f1", "kappa", "accuracy")
# métrica empleada para calcular el fitness.
#
# metodo_seleccion: ("ruleta", "rank", o "tournament").
# método para establecer la probabilidad de selección.
# Return
# ----------
# El índice ("integer") que ocupa el individuo seleccionado.
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------
if (!metodo_seleccion %in% c("ruleta", "rank", "tournament")) {
stop("El método de selección debe de ser ruleta, rank o tournament.")
}
# SELECCIÓN DE INDIVIDUOS
# ----------------------------------------------------------------------
# Se calcula la probabilidad de selección de cada individuo en función
# de su fitness.
if (metodo_seleccion == "ruleta") {
if (metrica %in% c("-mse", "-mae")) {
# Si el fitness es [-inf,0] se emplea 1/fitness
vector_fitness <- 1/vector_fitness
}
probabilidad_seleccion <- vector_fitness / sum(vector_fitness)
ind_seleccionado <- sample(
x = 1:length(vector_fitness),
size = 1,
prob = probabilidad_seleccion,
replace = TRUE
)
}else if (metodo_seleccion == "rank") {
# La probabilidad con este método es inversamente proporcional a la
# posición en la que quedan ordenados los individuos de menor a mayor
# fitness.
if (metrica %in% c("-mse", "-mae")) {
vector_fitness <- (-1/vector_fitness)
}
probabilidad_seleccion <- 1 / rank(-1*vector_fitness)
probabilidad_seleccion <- probabilidad_seleccion / sum(probabilidad_seleccion)
ind_seleccionado <- sample(
x = 1:length(vector_fitness),
size = 1,
prob = probabilidad_seleccion
)
}else if (metodo_seleccion == "tournament") {
# Se seleccionan aleatoriamente dos parejas de individuos.
ind_candidatos_a <- sample(x = 1:length(vector_fitness), size = 2)
ind_candidatos_b <- sample(x = 1:length(vector_fitness), size = 2)
# De cada pareja se selecciona el de mayor fitness.
ind_ganador_a <- ifelse(
vector_fitness[ind_candidatos_a[1]] > vector_fitness[ind_candidatos_a[2]],
ind_candidatos_a[1],
ind_candidatos_a[2]
)
ind_ganador_b <- ifelse(
vector_fitness[ind_candidatos_b[1]] > vector_fitness[ind_candidatos_b[2]],
ind_candidatos_b[1],
ind_candidatos_b[2]
)
# Se comparan los dos ganadores de cada pareja.
ind_seleccionado <- ifelse(
vector_fitness[ind_ganador_a] > vector_fitness[ind_ganador_b],
ind_ganador_a,
ind_ganador_b
)
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("----------------------", "\n")
cat("Individuo seleccionado", "\n")
cat("----------------------", "\n")
cat("Método selección:", metodo_seleccion, "\n")
cat("Índice seleccionado:", ind_seleccionado, "\n")
}
return(ind_seleccionado)
}Ejemplo
Se compara el patrón de selección del mejor individuo entre los métodos ruleta, rank y tournament. En primer lugar, se muestra un caso en el que la diferencia entre el mayor y el menor de los fitness no es muy acusada, y un segundo caso en el que sí lo es. En este ejemplo, el fitness se corresponde con el valor negativo del MSE, por lo que, cuanto más próximo a cero, mayor debe de ser la probabilidad de selección.
library(tidyverse)
fitness_poblacion <- c(-20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10)
selecciones_ruleta <- rep(NA, times=1000)
for (i in 1:1000) {
selecciones_ruleta[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "ruleta",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_ruleta <- data.frame(seleccion = selecciones_ruleta) %>%
mutate(metodo_seleccion = "ruleta")
selecciones_rank <- rep(NA, times=1000)
for (i in 1:1000) {
selecciones_rank[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "rank",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_rank <- data.frame(seleccion = selecciones_rank) %>%
mutate(metodo_seleccion = "rank")
selecciones_tournament <- rep(NA, times = 1000)
for (i in 1:1000) {
selecciones_tournament[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "tournament",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_tournament <- data.frame(seleccion = selecciones_tournament) %>%
mutate(metodo_seleccion = "tournament")
bind_rows(selecciones_ruleta, selecciones_rank, selecciones_tournament) %>%
ggplot(aes(x = seleccion)) +
geom_bar(fill = "#3182bd") +
scale_x_continuous(breaks = 1:11) +
labs(x="índice seleccionado") +
facet_grid(.~metodo_seleccion) +
theme_bw()
Cuando no existe una gran diferencia entre el individuo de mayor fitness y el resto, con el método rank, el individuo con mayor fitness se selecciona con mucha más frecuencia que el resto. Con los otros dos métodos, la probabilidad de selección decae de forma más gradual.
fitness_poblacion <- c(-200, -19, -18, -17, -16, -15, -14, -13, -12, -11, -1)
selecciones_ruleta <- rep(NA, times=1000)
for (i in 1:1000) {
selecciones_ruleta[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "ruleta",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_ruleta <- data.frame(seleccion = selecciones_ruleta) %>%
mutate(metodo_seleccion = "ruleta")
selecciones_rank <- rep(NA, times=1000)
for (i in 1:1000) {
selecciones_rank[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "rank",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_rank <- data.frame(seleccion = selecciones_rank) %>%
mutate(metodo_seleccion = "rank")
selecciones_tournament <- rep(NA, times = 1000)
for (i in 1:1000) {
selecciones_tournament[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "tournament",
metrica = "-mse",
verbose = FALSE
)
}
selecciones_tournament <- data.frame(seleccion = selecciones_tournament) %>%
mutate(metodo_seleccion = "tournament")
bind_rows(selecciones_ruleta, selecciones_rank, selecciones_tournament) %>%
ggplot(aes(x = seleccion)) +
geom_bar(fill = "#3182bd") +
scale_x_continuous(breaks = 1:11) +
labs(x="índice seleccionado") +
facet_grid(.~metodo_seleccion) +
theme_bw()
Cuando existe una gran diferencia entre el individuo de mayor fitness y el resto (uno o varios órdenes de magnitud), con el método ruleta, el individuo con mayor fitness se selecciona con mucha más frecuencia que el resto. A diferencia del caso anterior, en esta situación, la probabilidad de selección decae de forma más gradual con los métodos rank y tournament.
Teniendo en cuenta los comportamientos de selección de cada método, el método tournament parece ser la opción más equilibrada.
Cruzar dos individuos (crossover, recombinación)
El objetivo de esta etapa es generar, a partir de individuos ya existentes (parentales), nuevos individuos (descendencia) que combinen las características de los anteriores. Este es otro de los puntos del algoritmo en los que se puede seguir varias estrategias. Tres de las más empleadas son:
Cruzamiento a partir de uno solo punto: se selecciona aleatoriamente una posición que actúa como punto de corte. Cada individuo parental se divide en dos partes y se intercambian las mitades. Como resultado de este proceso, por cada cruce, se generan dos nuevos individuos.
Cruzamiento a partir múltiples puntos: se seleccionan aleatoriamente varias posiciones que actúan como puntos de corte. Cada individuo parental se divide por los puntos de corte y se intercambian las partes. Como resultado de este proceso, por cada cruce, se generan dos nuevos individuos.
Cruzamiento uniforme: el valor que toma cada posición del nuevo individuo se obtiene de uno de los dos parentales. Por lo general, la probabilidad de que el valor proceda de cada parental es la misma, aunque podría, por ejemplo, estar condicionada al fitness de cada uno. A diferencia de las anteriores estrategias, con esta, de cada cruce se genera un único descendiente.
cruzar_individuos <- function(parental_1,
parental_2,
metodo_cruce = "uniforme",
verbose = TRUE) {
# Esta función devuelve un individuo resultado de cruzar dos individuos
# parentales con el método de cruzamiento uniforme.
#
# Parameters
# ----------
# parental_1:
# vector que representa a un individuo.
# parental_2:
# vector que representa a un individuo.
# metodo_cruce: (uniforme", "punto_simple")
# estrategia de cruzamiento.
# Return
# ----------
# Un vector que representa a un nuevo individuo.
# Para crear el nuevo individuo, se emplea el método de cruzamiento uniforme,
# con la misma probabilidad de que el valor proceda de cada parental.
if (length(parental_1) != length(parental_2)) {
stop(paste0(
"La longitud de los dos vectores que representan a los ",
"individuos debe ser la misma."
))
}
if (!(metodo_cruce %in% c("uniforme", "punto_simple"))) {
stop("El método de cruzamiento debe ser: uniforme o punto_simple.")
}
# Se crea el vector que representa el nuevo individuo
descendencia <- rep(NA, times = length(parental_1))
if (metodo_cruce == "uniforme") {
# Se seleccionan aleatoriamente las posiciones que se heredan del parental_1.
herencia_parent_1 <- sample(
x = c(TRUE, FALSE),
size = length(parental_1),
replace = TRUE
)
# El resto de posiciones se heredan del parental_2.
herencia_parent_2 <- !herencia_parent_1
descendencia[herencia_parent_1] <- parental_1[herencia_parent_1]
descendencia[herencia_parent_2] <- parental_2[herencia_parent_2]
} else {
punto_cruce <- sample(
x = 2:length(parental_1),
size = 1
)
descendencia <- c(
parental_1[1:(punto_cruce - 1)],
parental_2[punto_cruce:length(parental_1)]
)
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("---------------", "\n")
cat("Cruce realizado", "\n")
cat("---------------", "\n")
cat("Método:", metodo_cruce, "\n")
cat("Parental 1", "\n")
cat(parental_1, "\n")
cat("Parental 2", "\n")
cat(parental_2, "\n")
cat("Descendencia", "\n")
cat(descendencia, "\n")
cat("\n")
}
return(descendencia)
}Ejemplo
Se obtiene un nuevo individuo a partir del cruce de los individuos c(T, T, T, T, T) y c(F, F, F, F, F, F).
cruzar_individuos(
parental_1 = c(T, T, T, T, T),
parental_2 = c(F, F, F, F, F),
metodo = "uniforme",
verbose = TRUE
)## ---------------
## Cruce realizado
## ---------------
## Método: uniforme
## Parental 1
## TRUE TRUE TRUE TRUE TRUE
## Parental 2
## FALSE FALSE FALSE FALSE FALSE
## Descendencia
## FALSE TRUE FALSE TRUE FALSE## [1] FALSE TRUE FALSE TRUE FALSEMutar individuo
Tras generar cada nuevo individuo de la descendencia, este se somete a un proceso de mutación en el que, cada una de sus posiciones, puede verse modificada con una probabilidad p. Este paso es importante para añadir diversidad al proceso y evitar que el algoritmo caiga en mínimos locales por que todos los individuos sean demasiado parecidos de una generación a otra.
mutar_individuo <- function(individuo, prob_mut = 0.01) {
# Este método somete al individuo a un proceso de mutación en el que, cada
# una de sus posiciones, puede verse modificada con una probabilidad `prob_mut`.
#
# Parameters
# ----------
# individuo:
# vector que representa a un individuo.
#
# prob_mut:
# probabilidad que tiene cada posición del vector de mutar.
#
# Return
# ----------
# Un vector que representa al individuo tras someterse a las mutaciones.
# Selección de posiciones a mutar.
posiciones_mutadas <- runif(n = length(individuo), min = 0, max = 1) < prob_mut
# Se modifica el valor de aquellas posiciones que hayan sido seleccionadas para
# mutar. Si el valor de prob_mut es muy bajo, las mutaciones serán muy poco
# frecuentes y el individuo devuelto será casi siempre igual al original.
individuo[posiciones_mutadas] <- !(individuo[posiciones_mutadas])
return(individuo)
}Ejemplo
Se somete a un individuo al proceso de mutación, con una probabilidad de mutación de 0.2.
mutar_individuo(individuo = c(T,T,T,T,T,T,T,T,T,T), prob_mut = 0.2)## [1] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSEComo resultado del cruce y mutación, el individuo generado puede exceder el número de predictores incluidos (n_max). Si se quiere evitar, se puede modificar la secuencia del individuo para que contenga como máximo n_max valores TRUE.
forzar_n_max = function(individuo, n_max){
# Este método modifica la secuencia del individuo de forma que, como máximo
# contenga n_max valores TRUE.
#
# Parameters
# ----------
# individuo:
# vector que representa a un individuo.
#
# n_max:
# número máximo de posiciones TRUE que puede tener el individuo.
# Return
# ----------
# Un vector que representa al individuo tras someterse a la corrección
# Se identifica si el número de TRUE es la secuencia supera a n_max.
n_exceso <- sum(individuo) - n_max
if (n_exceso > 0) {
# Se seleccionan aleatoriamente n_max posiciones con valor TRUE en
# la secuencia del individuo.
indices <- sample(
x = which(individuo == TRUE),
size = n_exceso,
replace = FALSE
)
individuo[indices] <- !(individuo[indices])
}
return(individuo)
}Algoritmo completo
En cada uno de los apartados anteriores se ha definido una de las etapas del algoritmo genético. A continuación, se combinan todas ellas dentro de una única función.
seleccionar_predictores_ga <- function(
x,
y,
n_poblacion = 20,
n_generaciones = 10,
n_max_predictores = NULL,
n_max_estricto = FALSE,
n_min_predictores = NULL,
modelo = "lm",
cv = 5,
metrica = NULL,
nivel_referencia = NULL,
elitismo = 0.1,
prob_mut = 0.01,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
n_tree = 50,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = FALSE,
verbose = TRUE,
...) {
# Parameters
# ----------
# x: "matrix" o "data.frame"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# n_poblacion: "integer"
# número total de individuos de la población.
#
# n_variables: "integer"
# longitud de los individuos.
#
# n_max_predictores: "integer"
# número máximo de predictores que puede contener un individuo.
#
# n_min_predictores: "integer"
# número mínimo de predictores que puede contener un individuo.
#
# modelo: ("lm", "glm", "randomforest").
# tipo de modelo empleado para calcular el fitness.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# metrica: ("-mse", "-mae", "f1", "kappa", "accuracy")
# métrica utilizada para calcular el fitness. Por defecto es el -mse para
# regresión y accuracy para clasificación. En el caso de clasificación
# binaria, se acepta también f1 y kappa.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# elitismo: "numeric"
# porcentaje de mejores individuos de la población actual que pasan
# directamente a la siguiente población. De esta forma, se asegura que, la
# siguiente generación, no sea nunca peor.
#
# prob_mut: "numeric"
# probabilidad que tiene cada posición del individuo de mutar.
#
# metodo_seleccion: ("ruleta", "rank", "tournament")
# método de selección de los individuos para los cruces.
#
# metodo_cruce: (uniforme", "punto_simple")
# estrategia de cruzamiento.
#
# parada_temprana: "logical"
# si durante las últimas `rondas_parada` generaciones el porcentaje de mejora
# entre mejores individuos no es superior al valor de `tolerancia_parada`,
# se detiene el algoritmo y no se crean nuevas generaciones.
#
# rondas_parada: "int"
# número de generaciones consecutivas sin mejora mínima para que se active
# la parada temprana.
#
# tolerancia_parada: "numeric"
# porcentaje de mejora que deben tener las generaciones consecutivas para
# considerar que hay cambio.
#
# n_tree: "integer"
# número de árboles del modelo randomforest.
#
# paralelizado: "logical"
# si se paraleliza o no el algoritmo de busqueda.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# Objeto "lista" con los resultados del algoritmo genético:
# fitness: fitness del mejor individuo de cada generación.
# mejor_individuo: información del mejor individuo de cada generación.
# diferencia_abs: mejora respecto a la generación anterior.
# porcentaje_mejora: porcentaje de mejora entre cada generación.
# df_resultados: dataframe con los resultados del proceso en cada generación.
# mejor_individuo_global: mejor individuo encontrado en todo el proceso.
library(foreach)
library(doParallel)
library(randomForest)
library(ranger)
library(caret)
library(MLmetrics)
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los parameters
# rondas_parada y tolerancia_parada.
if (parada_temprana && (is.null(rondas_parada) || is.null(tolerancia_parada))) {
stop(paste(
"Para activar la parada temprana es necesario indicar un valor",
"de rondas_parada y de tolerancia_parada."
))
}
start_time <- Sys.time()
# ALMACENAMIENTO DE RESULTADOS
# ----------------------------------------------------------------------------
# Por cada generación se almacena el mejor individuo, su fitness, y el porcentaje
# de mejora respecto a la última generación y la diferencia absoluta de mejora.
resultados_fitness <- vector(mode = "list", length = n_generaciones)
resultados_individuo <- vector(mode = "list", length = n_generaciones)
diferencia_abs <- vector(mode = "list", length = n_generaciones)
porcentaje_mejora <- vector(mode = "list", length = n_generaciones)
# CREACIÓN DE LA POBLACIÓN INICIAL
# ----------------------------------------------------------------------------
poblacion <- crear_poblacion(
n_poblacion = n_poblacion,
n_variables = ncol(x),
n_max = n_max_predictores,
n_min = n_min_predictores,
verbose = verbose
)
# ITERACIÓN DE POBLACIONES
# ----------------------------------------------------------------------------
for (i in 1:n_generaciones) {
if (verbose) {
cat("----------------", "\n")
cat("Generación:", i, "\n")
cat("----------------", "\n")
}
# CALCULAR FITNESS DE LOS INDIVIDUOS DE LA POBLACIÓN
# --------------------------------------------------------------------------
if (!paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
x = x,
y = y,
modelo = modelo,
cv = cv,
metrica = metrica,
nivel_referencia = nivel_referencia,
verbose = verbose
)
}
if (paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion_paral(
poblacion = poblacion,
x = x,
y = y,
modelo = modelo,
cv = cv,
metrica = metrica,
nivel_referencia = nivel_referencia,
verbose = verbose
)
}
# SE ALMACENA EL MEJOR INDIVIDUO DE LA POBLACIÓN ACTUAL
# --------------------------------------------------------------------------
fitness_mejor_individuo <- max(fitness_ind_poblacion)
mejor_individuo <- poblacion[which.max(fitness_ind_poblacion), ]
resultados_fitness[[i]] <- fitness_mejor_individuo
resultados_individuo[[i]] <- colnames(x)[mejor_individuo]
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA GENERACIÓN ANTERIOR
# ==========================================================================
# La diferencia solo puede calcularse a partir de la segunda generación.
if (i > 1) {
diferencia_abs[[i]] <- abs(resultados_fitness[[i]]- resultados_fitness[[i-1]])
}
# SE CALCULA LA MEJORA RESPECTO A LA GENERACIÓN ANTERIOR
# ==========================================================================
# La mejora solo puede calcularse a partir de la segunda generación.
if (i > 1) {
if(metrica %in% c("f1", "kappa", "accuracy")){
# Si el fitness se corresponde con una métrica con escala [0, +inf)
porcentaje_mejora[[i]] <- (resultados_fitness[[i]]-resultados_fitness[[i-1]]) /
resultados_fitness[[i-1]]
}
if(metrica %in% c("-mse", "-mae")){
# Si el fitness se corresponde con una métrica con escala (-inf, 0]
porcentaje_mejora[[i]] <- (resultados_fitness[[i-1]]-resultados_fitness[[i]]) /
resultados_fitness[[i]]
}
}
# NUEVA POBLACIÓN
# ==========================================================================
nueva_poblacion <- matrix(
data = NA,
nrow = nrow(poblacion),
ncol = ncol(poblacion)
)
# ELITISMO
# ==========================================================================
# El elitismo indica el porcentaje de mejores individuos de la población
# actual que pasan directamente a la siguiente población. De esta forma, se
# asegura que, la siguiente generación, no sea nunca inferior.
if (elitismo > 0) {
n_elitismo <- ceiling(nrow(poblacion) * elitismo)
posicion_n_mejores <- order(fitness_ind_poblacion, decreasing = TRUE)
posicion_n_mejores <- posicion_n_mejores[1:n_elitismo]
nueva_poblacion[1:n_elitismo, ] <- poblacion[posicion_n_mejores, ]
} else {
n_elitismo <- 0
}
# CREACIÓN DE NUEVOS INDIVIDUOS POR CRUCES
# ==========================================================================
for (j in (n_elitismo + 1):nrow(nueva_poblacion)) {
# Seleccionar parentales
indice_parental_1 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metrica = metrica,
metodo_seleccion = metodo_seleccion,
verbose = verbose
)
indice_parental_2 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metrica = metrica,
metodo_seleccion = metodo_seleccion,
verbose = verbose
)
parental_1 <- poblacion[indice_parental_1, ]
parental_2 <- poblacion[indice_parental_2, ]
# Cruzar parentales para obtener la descendencia
descendencia <- cruzar_individuos(
parental_1 = parental_1,
parental_2 = parental_2,
metodo_cruce = metodo_cruce,
verbose = verbose
)
# Mutar la descendencia
descendencia <- mutar_individuo(
individuo = descendencia,
prob_mut = prob_mut
)
# Si n_max_estricto=TRUE, se elimina el exceso de TRUEs en la
# secuencia de la descendencia.
if (n_max_estricto) {
descendencia <- forzar_n_max(
individuo = descendencia,
n_max = n_max_predictores
)
}
# Almacenar el nuevo descendiente en la nueva población: puede ocurrir que
# el individuo resultante del cruce y posterior mutación no contenga ningún
# predictor (todas sus posiciones son FALSE). Si esto ocurre, y para evitar
# que se produzca un error al ajustar el modelo, se introduce aleatoriamente
# un valor TRUE.
if (all(descendencia == FALSE)) {
descendencia[sample(x = 1:length(descendencia), size = 1)] <- TRUE
}
nueva_poblacion[j, ] <- descendencia
}
poblacion <- nueva_poblacion
# CRITERIO DE PARADA
# ==========================================================================
# Si durante las últimas n generaciones el mejor individuo no ha mejorado por
# una cantidad superior a tolerancia_parada, se detiene el algoritmo y no se
# crean nuevas generaciones.
if (parada_temprana && (i > rondas_parada)) {
ultimos_n <- tail(unlist(porcentaje_mejora), n = rondas_parada)
if (all(ultimos_n < tolerancia_parada)) {
print(paste(
"Algoritmo detenido en la generacion", i,
"por falta de mejora mínima de", tolerancia_parada,
"durante", rondas_parada,
"generaciones consecutivas."
))
break()
}
}
}
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------------
indice_mejor_individuo_global <- which.max(unlist(resultados_fitness))
mejor_fitness_global <- resultados_fitness[[indice_mejor_individuo_global]]
mejor_individuo_global <- resultados_individuo[[indice_mejor_individuo_global]]
# RESULTADOS
# ----------------------------------------------------------------------------
# La función devuelve una lista con 4 elementos:
# fitness:
# un vector con el fitness del mejor individuo de cada generación.
# mejor_individuo:
# una lista con la combinación de predictores del mejor individuo de cada
# generación.
# diferencia_abs:
# un vector con la diferencia absoluta entre el mejor individuo de cada
# generación.
# porcentaje_mejora:
# un vector con el porcentaje de mejora entre el mejor individuo de cada
# generación.
# df_resultados:
# un dataframe con todos los resultados anteriores.
# mejor_individuo_global:
# la combinación de predictores del mejor individuo encontrado en todo el
# proceso.
resultados_fitness <- unlist(resultados_fitness)
diferencia_abs <- c(NA, unlist(diferencia_abs))
porcentaje_mejora <- c(NA, unlist(porcentaje_mejora))
# Para poder añadir al dataframe la secuencia de predictores, se concatenan.
predictores <- resultados_individuo[!sapply(resultados_individuo, is.null)]
predictores <- sapply(
X = predictores,
FUN = function(x) {
paste(x, collapse = ", ")
}
)
df_resultados <- data.frame(
generacion = seq_along(resultados_fitness),
fitness = resultados_fitness,
predictores = predictores,
porcentaje_mejora = porcentaje_mejora,
diferencia_abs = diferencia_abs
)
# Frecuencia relativa con la que aparece cada predictor en el mejor individuo
# a lo largo de las generaciones.
frecuencia_seleccion <- tibble::enframe(
x = table(unlist(resultados_individuo)),
name = "pred",
value = "frec"
)
frecuencia_seleccion <- frecuencia_seleccion[order(-frecuencia_seleccion$frec), ]
frecuencia_seleccion$frec <- as.numeric(100*(frecuencia_seleccion$frec/i))
top_5_predictores <- head(frecuencia_seleccion, 5)$pred
resultados <- list(
fitness = resultados_fitness,
mejor_individuo = resultados_individuo,
diferencia_abs = diferencia_abs,
porcentaje_mejora = porcentaje_mejora,
df_resultados = df_resultados,
mejor_individuo_global = mejor_individuo_global,
frecuencia_seleccion = frecuencia_seleccion,
top_5_predictores = top_5_predictores
)
end_time <- Sys.time()
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(resultados, "class") <- "seleccion_ga"
attr(resultados, 'fecha_creacion') <- end_time
attr(resultados, 'duracion_seleccion') <- difftime(end_time, start_time)
attr(resultados, 'generaciones') <- i
attr(resultados, 'predictores_seleccionados') <- mejor_individuo_global
attr(resultados, 'top_5_predictores') <- top_5_predictores
attr(resultados, 'fitness_optimo') <- mejor_fitness_global
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'modelo') <- modelo
attr(resultados, 'metrica') <- metrica
attr(resultados, 'elitismo') <- elitismo
attr(resultados, 'prob_mut') <- prob_mut
attr(resultados, 'metodo_seleccion') <- metodo_seleccion
attr(resultados, 'metodo_cruce') <- metodo_cruce
attr(resultados, 'parada_temprana') <- parada_temprana
attr(resultados, 'rondas_parada') <- rondas_parada
attr(resultados, 'tolerancia_parada') <- tolerancia_parada
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
cat("--------------------", "\n")
cat("Selección finalizada", "\n")
cat("--------------------", "\n")
cat("Fecha finalización:", as.character(Sys.time()), "\n")
cat("Duración selección: ")
print(difftime(end_time, start_time))
cat("Número de generaciones:", i, "\n")
cat("Predictores seleccionados:", mejor_individuo_global,"\n")
cat("Fitness final:", mejor_fitness_global, "\n")
cat("Top 5 predictores:", top_5_predictores,"\n")
cat("Modelo empleado:", modelo, "\n")
cat("Métrica empleada:", metrica, "\n")
cat("\n")
return(resultados)
}Se añade también un método print para la clase seleccion_ga.
print.seleccion_ga <- function(obj){
cat("-------------------------------", "\n")
cat("Selección de predictores por GA", "\n")
cat("-------------------------------", "\n")
cat("Fecha creación:", as.character(attr(obj, 'fecha_creacion')), "\n")
cat("Predictores seleccionados:", attr(obj, 'predictores_seleccionados'),"\n")
cat("Fitness final:", attr(obj, 'fitness_optimo'), "\n")
cat("Top 5 predictores:", attr(obj, 'top_5_predictores'),"\n")
cat("Duración selección:", attr(obj, 'duracion_seleccion'),
attr(attr(obj, 'duracion_seleccion'), "units"),"\n")
cat("Número de generaciones:", attr(obj, 'generaciones'), "\n")
cat("Tamaño población:", attr(obj, 'n_poblacion'), "\n")
cat("Modelo empleado:", attr(obj, 'modelo'), "\n")
cat("Métrica empleada:", attr(obj, 'metrica'), "\n")
cat("Probabilidad mutación:", attr(obj, 'prob_mut'), "\n")
cat("Elitismo:", attr(obj, 'elitismo'), "\n")
cat("Método selección:", attr(obj, 'metodo_seleccion'), "\n")
cat("Método cruce:", attr(obj, 'metodo_cruce'), "\n")
cat("\n")
}Ejemplo regresión
Simulación de datos
En el paquete mlbench puede encontrarse, entre muchos otros data sets, el problema de regresión propuesto por Friedman (1991) y Breiman (1996). Con la función mlbench.friedman1() se puede generar un conjunto de datos simulados que siguen la ecuación:
\[y = 10 sin(\pi x_1 x_2) + 20(x_3 − 0.5)^2 + 10x_4 + 5x5 + \epsilon\]
Además de las primeras 5 columnas (\(x_1,...,x_5\)), que están relacionadas con la variable respuesta, se añaden automáticamente 5 columnas adicionales (\(x_6,...,x_10\)) que siguen una distribución uniforme [0,1] y que, por lo tanto, no guardan relación alguna con la variable respuesta.
library(mlbench)
library(dplyr)
set.seed(123)
simulacion <- mlbench.friedman1(n = 500, sd = 1)
# El objeto simulación es una lista que contiene una matriz con los 10 predictores y
# un vector con la variable respuesta. Se unen todos los datos en un único dataframe.
datos <- as.data.frame(simulacion$x)
datos$y <- simulacion$y
colnames(datos) <- c(paste("x", 1:10, sep = ""), "y")
datos <- datos %>% select(y, everything())
head(datos)Se añaden además 20 columnas adicionales con valores aleatorios distribuidos de forma normal.
n <- 500
p <- 20
ruido <- matrix(rnorm(n * p), nrow = n)
colnames(ruido) <- paste("x", 11:(10+p), sep = "")
ruido <- as.data.frame(ruido)
datos <- bind_cols(datos, ruido)
head(datos)Como resultado se ha generado un conjunto de datos con 30 predictores, de los cuales, solo 5 están realmente relacionados con la variable respuesta.
Selección de predictores
En primer lugar se emplean pocos individuos y generaciones para mostrar la información que se imprime por pantalla cuando verbose = TRUE.
seleccion <- seleccionar_predictores_ga(
x = datos[, -1],
y = datos$y,
n_poblacion = 3,
n_generaciones = 2,
n_max_predictores = 5,
n_max_estricto = FALSE,
n_min_predictores = 1,
elitismo = 0.01,
prob_mut = 0.1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
modelo = "randomforest",
n_tree = 50,
metrica = "-mse",
cv = 5,
parada_temprana = TRUE,
rondas_parada = 5,
tolerancia_parada = 0.01,
verbose = TRUE
)## Población inicial creada
## ------------------------
## Fecha creación: 2020-02-01 18:58:14
## Número de individuos = 3
## Número de predictores mínimo por individuo = 1
## Número de predictores máximo por individuo = 5
##
## ----------------
## Generación: 1
## ----------------
## Individuo 1 : FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -18.50789
## Fitness: -18.50789
## Modelo de evaluación: Random Forest
##
## Individuo 2 : FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -20.12239
## Fitness: -20.12239
## Modelo de evaluación: Random Forest
##
## Individuo 3 : FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -31.47203
## Fitness: -31.47203
## Modelo de evaluación: Random Forest
##
## ------------------
## Población evaluada
## ------------------
## Modelo empleado para el cálculo del fitness: randomforest
## Métrica empleado para el cálculo del fitness: -mse
## Mejor fitness encontrado: -18.50789
## Mejor secuencia encontrada:
## [1] FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE TRUE FALSE FALSE
## Mejores predictores encontrados:
## [1] "x4" "x8" "x28"
##
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 2
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 1
## ---------------
## Cruce realizado
## ---------------
## Método: uniforme
## Parental 1
## FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Parental 2
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## Descendencia
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 1
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 2
## ---------------
## Cruce realizado
## ---------------
## Método: uniforme
## Parental 1
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## Parental 2
## FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Descendencia
## FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
##
## ----------------
## Generación: 2
## ----------------
## Individuo 1 : FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -18.50789
## Fitness: -18.50789
## Modelo de evaluación: Random Forest
##
## Individuo 2 : FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -18.95568
## Fitness: -18.95568
## Modelo de evaluación: Random Forest
##
## Individuo 3 : FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## El individuo ha sido evaluado
## -----------------------------
## Métrica: -mse
## Valor métrica: -25.43676
## Fitness: -25.43676
## Modelo de evaluación: Random Forest
##
## ------------------
## Población evaluada
## ------------------
## Modelo empleado para el cálculo del fitness: randomforest
## Métrica empleado para el cálculo del fitness: -mse
## Mejor fitness encontrado: -18.50789
## Mejor secuencia encontrada:
## [1] FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE TRUE FALSE FALSE
## Mejores predictores encontrados:
## [1] "x4" "x8" "x28"
##
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 2
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 1
## ---------------
## Cruce realizado
## ---------------
## Método: uniforme
## Parental 1
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Parental 2
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## Descendencia
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 1
## ----------------------
## Individuo seleccionado
## ----------------------
## Método selección: tournament
## Índice seleccionado: 2
## ---------------
## Cruce realizado
## ---------------
## Método: uniforme
## Parental 1
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## Parental 2
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Descendencia
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
##
## --------------------
## Selección finalizada
## --------------------
## Fecha finalización: 2020-02-01 18:58:15
## Duración selección: Time difference of 0.4845126 secs
## Número de generaciones: 2
## Predictores seleccionados: x4 x8 x28
## Fitness final: -18.50789
## Top 5 predictores: x28 x4 x8
## Modelo empleado: randomforest
## Métrica empleada: -mseseleccion## -------------------------------
## Selección de predictores por GA
## -------------------------------
## Fecha creación: 2020-02-01 18:58:15
## Predictores seleccionados: x4 x8 x28
## Fitness final: -18.50789
## Top 5 predictores: x28 x4 x8
## Duración selección: 0.4845126 secs
## Número de generaciones: 2
## Tamaño población: 3
## Modelo empleado: randomforest
## Métrica empleada: -mse
## Probabilidad mutación: 0.1
## Elitismo: 0.01
## Método selección: tournament
## Método cruce: uniforme# Mejor fitness de cada generación.
seleccion$fitness## [1] -18.50789 -18.50789# Mejor individuo de cada generación.
seleccion$mejor_individuo## [[1]]
## [1] "x4" "x8" "x28"
##
## [[2]]
## [1] "x4" "x8" "x28"# Differencia absoluta del mejor fitness entre generaciones.
seleccion$diferencia_abs## [1] NA 0# Resultados en forma de dataframe
seleccion$df_resultadosSe repite el proceso pero, esta vez, con 50 individuos por generación y 20 generaciones. Se especifica que verbose = FALSE para no saturar la consola.
seleccion <- seleccionar_predictores_ga(
x = datos[, -1],
y = datos$y,
n_poblacion = 50,
n_generaciones = 20,
n_max_predictores = 5,
n_max_estricto = FALSE,
n_min_predictores = 1,
elitismo = 0.01,
prob_mut = 0.1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
modelo = "randomforest",
n_tree = 50,
metrica = "-mse",
cv = 5,
parada_temprana = TRUE,
rondas_parada = 5,
tolerancia_parada = 0.01,
verbose = FALSE
)## [1] "Algoritmo detenido en la generacion 15 por falta de mejora mínima de 0.01 durante 5 generaciones consecutivas."
## --------------------
## Selección finalizada
## --------------------
## Fecha finalización: 2020-02-01 18:59:32
## Duración selección: Time difference of 1.286931 mins
## Número de generaciones: 15
## Predictores seleccionados: x1 x2 x3 x4 x5 x13 x23
## Fitness final: -5.076828
## Top 5 predictores: x4 x1 x2 x5 x3
## Modelo empleado: randomforest
## Métrica empleada: -mseseleccion## -------------------------------
## Selección de predictores por GA
## -------------------------------
## Fecha creación: 2020-02-01 18:59:32
## Predictores seleccionados: x1 x2 x3 x4 x5 x13 x23
## Fitness final: -5.076828
## Top 5 predictores: x4 x1 x2 x5 x3
## Duración selección: 1.286931 mins
## Número de generaciones: 15
## Tamaño población: 50
## Modelo empleado: randomforest
## Métrica empleada: -mse
## Probabilidad mutación: 0.1
## Elitismo: 0.01
## Método selección: tournament
## Método cruce: uniformeMejor individuo
El objeto devuelto por la función seleccionar_predictores almacena el mejor individuo de cada generación, la evolución del fitness a lo largo de las generaciones y el mejor individuo encontrado a lo largo de todo el proceso.
seleccion$mejor_individuo_global## [1] "x1" "x2" "x3" "x4" "x5" "x13" "x23"Evolución del error
La siguiente función extrae la información de las generaciones de un objeto de clase `` y crea un gráfico que muestra cómo evoluciona el fitness del mejor individuo a medida que avanzan las generaciones.
plot.seleccion_ga <- function(obj,
evolucion_fitness = TRUE,
frecuencia_seleccion = TRUE){
if (evolucion_fitness) {
p1 <- ggplot(data = obj$df_resultados, aes(x = generacion, y = fitness)) +
geom_line() +
geom_point() +
labs(title = 'Evolución del mejor individuo',
x = 'generacion',
y = 'fitness'
) +
theme_bw()
print(p1)
}
if (frecuencia_seleccion) {
p2 <- ggplot(data = obj$frecuencia_seleccion,
aes(x = reorder(pred, frec), y = frec, fill = frec)) +
geom_col() +
scale_fill_viridis_c() +
coord_flip() +
labs(title = 'Frecuencia de selección',
x = 'predictor',
y = 'frecuencia'
) +
theme_bw() +
theme(legend.position = "none")
print(p2)
}
}plot.seleccion_ga(
obj = seleccion,
evolucion_fitness = TRUE,
frecuencia_seleccion = FALSE
)
Frecuencia selección variables
Tal y como ocurre en este caso (sobre todo si se emplea ranfom forest con muchos árboles), el mejor individuo puede incorporar algunos predictores no relevantes o dejarse fuera alguno de los importantes. Una forma de ser críticos con el resultado obtenido al aplicar el algoritmo genético, es seguir la siguiente idea: si realmente existen variables útiles, es de esperar que estas tiendan a aparecer en el mejor individuo de cada generación, mientras que, las variables no relevantes, cambiarán con más frecuencia.
Como se han almacenado las variables que incluye el mejor individuo de cada generación, se puede calcular la frecuencia con la que ha aparecido cada una.
seleccion$frecuencia_seleccionplot.seleccion_ga(
obj = seleccion,
evolucion_fitness = FALSE,
frecuencia_seleccion = TRUE
)
Una forma de facilitar la valoración de estas frecuencias es compararlas frente a la frecuencia esperada solo por azar. La probabilidad de que una variable \(i\) perteneciente un conjunto de variables de tamaño \(M\) esté incluida en un subconjunto de tamaño \(N\) obtenido por muestreo aleatorio sin reposición es:
\[\text{P(i en N)} = N/M\]
Como el tamaño de cada individuo varía, hay que tenerlo en cuenta en el cálculo. Supóngase 5 individuos de tamaños \(K = \{5, 5, 3, 2, 1\}\) y que el total de variables disponibles es \(M\). El número de individuos en los que se espera que una determinada variable haya sido seleccionada por azar es:
\[\sum_{i \in K} \frac{i}{M}\]
La siguiente función realiza el cálculo de la frecuencia esperada.
calcular_frecuencia_esperada <- function(n_variables, individuos){
# Argumentos:
# n_variables: número total de variables disponibles.
# individuos: lista con el mejor individuo de cada generación.
# Retorno:
# Porcentaje de individuos en los que se espera que una variable aparezca
# solo por azar.
n_individuos <- sum(!sapply(X = individuos, FUN = is.null))
tamayos <- sapply(X = individuos, FUN = length)
tamayos <- tamayos[tamayos != 0]
frecuencia_esperada <- 0
for (i in tamayos) {
frecuencia_esperada <- frecuencia_esperada + i/n_variables
}
return(frecuencia_esperada/n_individuos)
}calcular_frecuencia_esperada(
n_variables = ncol(datos)-1,
individuos = seleccion$mejor_individuo
)## [1] 0.2866667Puede observarse que, las variables que aparecen con más frecuencia en los mejores individuos de cada generación, son las que están realmente relacionadas con la variable respuesta. Además, todas ellas superan la frecuencia esperada por azar. Combinando el resultado obtenido en el individuo final del algoritmo genético con la información de las frecuencias, se pueden obtener buenos resultados de selección.
Si se utiliza elitismo, el mejor individuo de una generación pasa por las siguientes n generaciones de forma directa. De ser así, las variables no relevantes que arrastre este individuo aparecerán con una frecuencia superior a la esperada por azar aun cuando no son predictores útiles.
Ejemplo clasificación
Datos
El set de datos titanic_train del paquete titanic contiene información sobre los pasajeros del transatlántico titanic así como si sobrevivieron o no al naufragio. Se pretende identificar las variables con las que mejor se puede predecir si un pasajero vivió o falleció.
library(titanic)
# Se eliminan para este ejemplo 3 variables.
datos <- titanic_train %>% select(-Name, -Ticket, -Cabin)
# Se eliminan observaciones incompletas.
datos <- na.omit(datos)
# Se convierten a factor las variables de tipo character.
datos <- datos %>% mutate_if(is.character, as.factor)
datos <- datos %>% mutate(Survived = as.factor(Survived))Selección de predictores
seleccion <- seleccionar_predictores_ga(
x = datos[, -2],
y = datos[, 2],
n_poblacion = 50,
n_generaciones = 20,
n_max_predictores = 8,
n_min_predictores = 1,
elitismo = 0.01,
prob_mut = 0.01,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
modelo = "randomforest",
metrica = "f1",
nivel_referencia = "1",
cv = 5,
parada_temprana = TRUE,
rondas_parada = 5,
tolerancia_parada = 0.01,
verbose = FALSE
)## [1] "Algoritmo detenido en la generacion 6 por falta de mejora mínima de 0.01 durante 5 generaciones consecutivas."
## --------------------
## Selección finalizada
## --------------------
## Fecha finalización: 2020-02-01 18:59:57
## Duración selección: Time difference of 24.88772 secs
## Número de generaciones: 6
## Predictores seleccionados: Pclass Sex Age SibSp Fare
## Fitness final: 0.7794929
## Top 5 predictores: Age Fare Pclass Sex SibSp
## Modelo empleado: randomforest
## Métrica empleada: f1seleccion## -------------------------------
## Selección de predictores por GA
## -------------------------------
## Fecha creación: 2020-02-01 18:59:57
## Predictores seleccionados: Pclass Sex Age SibSp Fare
## Fitness final: 0.7794929
## Top 5 predictores: Age Fare Pclass Sex SibSp
## Duración selección: 24.88772 secs
## Número de generaciones: 6
## Tamaño población: 50
## Modelo empleado: randomforest
## Métrica empleada: f1
## Probabilidad mutación: 0.01
## Elitismo: 0.01
## Método selección: tournament
## Método cruce: uniformeEvolución del error
plot.seleccion_ga(
obj = seleccion,
evolucion_fitness = TRUE,
frecuencia_seleccion = FALSE
)
Paralelización
Uno de los inconvenientes de los algoritmos genéticos es su alto requerimiento computacional. Por ejemplo, si se establecen 20 generaciones con 50 individuos por generación y una estimación del fitness mediante validación cruzada de 5 particiones, en total, se tienen que ajustar \(50x20x5 = 5000\) modelos.
Dos de las estrategias que se pueden emplear para agilizar el proceso son:
Parada temprana: detener el algoritmo si tras n generaciones consecutivas no se ha conseguido mejora. Esta estrategia está implementada en los ejemplos anteriores.
Paralelización: ajustar varios modelos de forma simultánea empleando múltiples cores del ordenador.
El algoritmo genético propuesto en este documento tiene dos etapas que pueden ser paralelizadas: la evaluación de los individuos de una población (calcular_fitness_poblacion()) y la validación cruzada de cada individuo (calcular_fitness_individuo()). Se procede a paralelizar la primera y comparar el impacto en el tiempo de ejecución.
Versión paralelizada
ESTA IMPLEMENTACIÓN NO FUNCIONA EN WINDOWS
calcular_fitness_poblacion_paral <- function(poblacion,
x,
y,
modelo,
cv,
seed = 123,
metrica = NULL,
nivel_referencia = NULL,
verbose = TRUE) {
# Parameters
# ----------
# poblacion: "matrix"
# matriz boleana que representa la población de individuos.
#
# x: "matrix" o "data.frame"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# modelo: ("lm", "glm", "randomforest").
# tipo de modelo empleado para calcular el fitness.
#
# metrica: ("-mse", "-mae", "f1", "kappa", "accuracy").
# métrica utilizada para calcular el fitness. Por defecto es el -mse para
# regresión y accuracy para clasificación. En el caso de clasificación
# binaria, se acepta también f1 y kappa.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
#
# Return
# ----------
# Vector ("numeric") con el fitness de todos los individuos de la población,
# obtenido por validación cruzada. El orden de los valores se corresponde con
# el orden de las filas de la matriz población.
# LIBRERIAS NECESARIAS
# ----------------------------------------------------------------------------
library(foreach)
library(doParallel)
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (!modelo %in% c("lm", "glm", "randomforest")) {
stop(
paste("El modelo empleado para calcular el fitness debe ser",
"lm (regresion lineal)",
"glm (regresión logistica)",
"o randomforest."
)
)
}
if (ncol(poblacion) != ncol(x)) {
stop(
"n_variables de la población debe ser igual al número de columnas de x."
)
}
# Si el modelo es de regresión lineal
if (modelo == "lm") {
if (!is.numeric(y)) {
stop(
"El método lm solo es valido para problemas de regresión."
)
}
if (!is.null(metrica)) {
if (!metrica %in% c("-mse", "-mae")) {
stop(
paste(
"La métrica", metrica, "no es válida para problemas de regresión.",
"Emplear el -mse o -mae."
)
)
}
}else{
# Por defecto se emplea -mse
metrica <- "-mse"
}
}
# Si el modelo es de regresión logística
if (modelo == "glm") {
if (!is.null(metrica)) {
if (!metrica %in% c("accuracy", "f1", "kappa")) {
stop(
"Las métricas permitidas con el modelo glm son: accuracy, f1 y kappa."
)
}
}else{
# Por defecto se emplea accuracy
metrica <- "accuracy"
}
if (!is.character(y) & !is.factor(y)) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación,",
"la variable respuesta debe ser character o factor."
)
)
}
if (is.character(y)) {
y <- as.factor(y)
}
if (length(levels(y)) != 2) {
stop(
paste(
"El modelo glm solo puede emplearse para problemas de clasificación binaria,",
"la variable respuesta debe tener dos niveles."
)
)
}
}
# Si el modelo es de random forest
if (modelo == "randomforest") {
if (!is.null(metrica)) {
if (metrica %in% c("-mse", "-mae") && !is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de regresión.")
)
}
if (metrica %in% c("accuracy", "f1", "kappa") && is.numeric(y)) {
stop(
paste(
"La métrica", metrica, "solo es válida para problemas de clasificación.")
)
}
if (metrica %in% c("f1", "kappa") && length(unique(y)) > 2) {
stop(
paste(
"La métrica",
metrica,
"solo es válida para problemas de clasificación binaria."
)
)
}
if (metrica %in% c("f1", "kappa") && is.null(nivel_referencia)) {
stop(
"Para la métrica kappa es necesario indicar el nivel de referencia."
)
}
}
if (is.null(metrica)) {
if (is.numeric(y)) {
metrica <- "-mse"
}else{
metrica <- "accuracy"
}
}
}
# SELECCIÓN FUNCIÓN AUXILIAR PARA CALCULAR EL FITNESS
# ----------------------------------------------------------------------------
# Tipo de modelo utilizado para calcular el fitness.
if (modelo == "lm") {
calcular_fitness_individuo <- calcular_fitness_individuo_lm
} else if (modelo == "randomforest") {
calcular_fitness_individuo <- calcular_fitness_individuo_rf
} else if (modelo == "glm") {
calcular_fitness_individuo <- calcular_fitness_individuo_glm
}
# CÁLCULO DEL FITNNES DE CADA INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
# Vector donde almacenar el fitness de cada individuo.
fitness_poblacion <- rep(NA, times = nrow(poblacion))
# Se emplean todos los cores del ordenador menos 1.
registerDoParallel(parallel::detectCores() - 1)
# Se emplea la función foreach para paralelizar.
fitness_poblacion <- foreach::foreach(i = 1:nrow(poblacion)) %dopar% {
individuo <- poblacion[i, ]
if (verbose) {
print(paste("Individuo", i, ":", paste(individuo, collapse = " ")))
}
fitness_individuo <- calcular_fitness_individuo(
x = x[, individuo, drop = FALSE],
y = y,
cv = cv,
seed = seed,
metrica = metrica,
nivel_referencia = nivel_referencia,
verbose = verbose
)
fitness_individuo
}
fitness_poblacion <- unlist(fitness_poblacion)
# Se vuelve a un único core.
options(cores = 1)
# MEJOR INDIVIDUO DE LA POBLACIÓN
# ------------------------------------------------------------------------
# Se identifica el mejor individuo de toda la población, el de mayor
# fitness.
indice_mejor_individuo <- which.max(fitness_poblacion)
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("------------------", "\n")
cat("Población evaluada", "\n")
cat("------------------", "\n")
cat("Modelo empleado para el cálculo del fitness:", modelo, "\n")
cat("Métrica empleado para el cálculo del fitness:", metrica, "\n")
cat("Mejor fitness encontrado:", fitness_poblacion[indice_mejor_individuo], "\n")
cat("Mejor secuencia encontrada:", "\n")
print(poblacion[indice_mejor_individuo,])
cat("Mejores predictores encontrados:", "\n")
print(colnames(x)[poblacion[indice_mejor_individuo,]])
cat("\n")
}
return(fitness_poblacion)
}Para incluir la opción de paralelizado, se repite la función del algoritmo completo, esta vez, incluyendo el argumento paralelizado con el que el usuario pueda especificar que se emplee la función calcular_fitness_poblacion o calcular_fitness_poblacion_paral.
seleccionar_predictores_ga <- function(
x,
y,
n_poblacion = 20,
n_generaciones = 10,
n_max_predictores = NULL,
n_max_estricto = FALSE,
n_min_predictores = NULL,
modelo = "lm",
cv = 5,
metrica = NULL,
nivel_referencia = NULL,
elitismo = 0.1,
prob_mut = 0.01,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
n_tree = 50,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = FALSE,
verbose = TRUE,
...) {
# Parameters
# ----------
# x: "matrix" o "data.frame"
# matriz de predictores.
#
# y: "vector"
# variable respuesta.
#
# n_poblacion: "integer"
# número total de individuos de la población.
#
# n_variables: "integer"
# longitud de los individuos.
#
# n_max_predictores: "integer"
# número máximo de predictores que puede contener un individuo.
#
# n_min_predictores: "integer"
# número mínimo de predictores que puede contener un individuo.
#
# modelo: ("lm", "glm", "randomforest").
# tipo de modelo empleado para calcular el fitness.
#
# cv: "integer"
# número de particiones de validación cruzada.
#
# seed: "integer"
# semilla para garantizar reproducibilidad en el proceso de CV.
#
# metrica: ("-mse", "mae", "f1", "kappa", "accuracy")
# métrica utilizada para calcular el fitness. Por defecto es el -mse para
# regresión y accuracy para clasificación. En el caso de clasificación
# binaria, se acepta también f1 y kappa.
#
# nivel_referencia: "character"
# valor de la variable respuesta considerado como referencia. Necesario
# cuando la métrica es f1 o kappa.
#
# elitismo: "numeric"
# porcentaje de mejores individuos de la población actual que pasan
# directamente a la siguiente población. De esta forma, se asegura que, la
# siguiente generación, no sea nunca peor.
#
# prob_mut: "numeric"
# probabilidad que tiene cada posición del individuo de mutar.
#
# metodo_seleccion: ("ruleta", "rank", "tournament")
# método de selección de los individuos para los cruces.
#
# metodo_cruce: (uniforme", "punto_simple")
# estrategia de cruzamiento.
#
# parada_temprana: "logical"
# si durante las últimas `rondas_parada` generaciones el porcentaje de mejora
# entre mejores individuos no es superior al valor de `tolerancia_parada`,
# se detiene el algoritmo y no se crean nuevas generaciones.
#
# rondas_parada: "int"
# número de generaciones consecutivas sin mejora mínima para que se active
# la parada temprana.
#
# tolerancia_parada: "numeric"
# porcentaje de mejora que deben las generaciones consecutivas para
# considerar que hay cambio.
#
# n_tree: "integer"
# número de árboles del modelo randomforest.
#
# paralelizado: "logical"
# si se paraleliza o no el algoritmo de busqueda.
#
# verbose: "logical"
# mostrar información del proceso por pantalla.
# Return
# ----------
# Objeto "lista" con los resultados del algoritmo genético:
# fitness: fitness del mejor individuo de cada generación.
# mejor_individuo: información del mejor individuo de cada generación.
# diferencia_abs: mejora respecto a la generación anterior.
# porcentaje_mejora: porcentaje de mejora entre cada generación.
# df_resultados: dataframe con los resultados del proceso en cada generación.
# mejor_individuo_global: mejor individuo encontrado en todo el proceso.
# LIBRERIAS NECESARIAS
# ----------------------------------------------------------------------------
library(foreach)
library(doParallel)
library(randomForest)
library(ranger)
library(caret)
library(MLmetrics)
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los parameters
# rondas_parada y tolerancia_parada.
if (parada_temprana && (is.null(rondas_parada) || is.null(tolerancia_parada))) {
stop(paste(
"Para activar la parada temprana es necesario indicar un valor",
"de rondas_parada y de tolerancia_parada."
))
}
if (paralelizado) {
if (Sys.info()["sysname"] == "Windows") {
cat("Windows no permite la ejecución en paralelo", "\n")
cat("Se ejecuta en modo secuencial", "\n")
paralelizado <- FALSE
}
}
start_time <- Sys.time()
# ALMACENAMIENTO DE RESULTADOS
# ----------------------------------------------------------------------------
# Por cada generación se almacena el mejor individuo, su fitness, y la mejora
# respecto a la última generación.
resultados_fitness <- vector(mode = "list", length = n_generaciones)
resultados_individuo <- vector(mode = "list", length = n_generaciones)
diferencia_abs <- vector(mode = "list", length = n_generaciones)
porcentaje_mejora <- vector(mode = "list", length = n_generaciones)
# CREACIÓN DE LA POBLACIÓN INICIAL
# ----------------------------------------------------------------------------
poblacion <- crear_poblacion(
n_poblacion = n_poblacion,
n_variables = ncol(x),
n_max = n_max_predictores,
n_min = n_min_predictores,
verbose = verbose
)
# ITERACIÓN DE POBLACIONES
# ----------------------------------------------------------------------------
for (i in 1:n_generaciones) {
if (verbose) {
cat("----------------", "\n")
cat("Generación:", i, "\n")
cat("----------------", "\n")
}
# CALCULAR FITNESS DE LOS INDIVIDUOS DE LA POBLACIÓN
# --------------------------------------------------------------------------
if (!paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
x = x,
y = y,
modelo = modelo,
cv = cv,
metrica = metrica,
nivel_referencia = nivel_referencia,
verbose = verbose
)
}
if (paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion_paral(
poblacion = poblacion,
x = x,
y = y,
modelo = modelo,
cv = cv,
metrica = metrica,
nivel_referencia = nivel_referencia,
verbose = verbose
)
}
# SE ALMACENA EL MEJOR INDIVIDUO DE LA POBLACIÓN ACTUAL
# --------------------------------------------------------------------------
fitness_mejor_individuo <- max(fitness_ind_poblacion)
mejor_individuo <- poblacion[which.max(fitness_ind_poblacion), ]
resultados_fitness[[i]] <- fitness_mejor_individuo
resultados_individuo[[i]] <- colnames(x)[mejor_individuo]
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA GENERACIÓN ANTERIOR
# ==========================================================================
# La diferencia solo puede calcularse a partir de la segunda generación.
if (i > 1) {
diferencia_abs[[i]] <- abs(resultados_fitness[[i]]- resultados_fitness[[i-1]])
}
# SE CALCULA LA MEJORA RESPECTO A LA GENERACIÓN ANTERIOR
# ==========================================================================
# La mejora solo puede calcularse a partir de la segunda generación.
if (i > 1) {
if(metrica %in% c("f1", "kappa", "accuracy")){
# Si el fitness se corresponde con una métrica con escala [0, +inf)
porcentaje_mejora[[i]] <- (resultados_fitness[[i]]-resultados_fitness[[i-1]]) /
resultados_fitness[[i-1]]
}
if(metrica %in% c("-mse", "-mae")){
# Si el fitness se corresponde con una métrica con escala (-inf, 0]
porcentaje_mejora[[i]] <- (resultados_fitness[[i-1]]-resultados_fitness[[i]]) /
resultados_fitness[[i]]
}
}
# NUEVA POBLACIÓN
# ==========================================================================
nueva_poblacion <- matrix(
data = NA,
nrow = nrow(poblacion),
ncol = ncol(poblacion)
)
# ELITISMO
# ==========================================================================
# El elitismo indica el porcentaje de mejores individuos de la población
# actual que pasan directamente a la siguiente población. De esta forma, se
# asegura que, la siguiente generación, no sea nunca inferior.
if (elitismo > 0) {
n_elitismo <- ceiling(nrow(poblacion) * elitismo)
posicion_n_mejores <- order(fitness_ind_poblacion, decreasing = TRUE)
posicion_n_mejores <- posicion_n_mejores[1:n_elitismo]
nueva_poblacion[1:n_elitismo, ] <- poblacion[posicion_n_mejores, ]
} else {
n_elitismo <- 0
}
# CREACIÓN DE NUEVOS INDIVIDUOS POR CRUCES
# ==========================================================================
for (j in (n_elitismo + 1):nrow(nueva_poblacion)) {
# Seleccionar parentales
indice_parental_1 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metrica = metrica,
metodo_seleccion = metodo_seleccion,
verbose = verbose
)
indice_parental_2 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metrica = metrica,
metodo_seleccion = metodo_seleccion,
verbose = verbose
)
parental_1 <- poblacion[indice_parental_1, ]
parental_2 <- poblacion[indice_parental_2, ]
# Cruzar parentales para obtener la descendencia
descendencia <- cruzar_individuos(
parental_1 = parental_1,
parental_2 = parental_2,
metodo_cruce = metodo_cruce,
verbose = verbose
)
# Mutar la descendencia
descendencia <- mutar_individuo(
individuo = descendencia,
prob_mut = prob_mut
)
# Si n_max_estricto=TRUE, se elimina el exceso de TRUEs en la
# secuencia de la descendencia.
if (n_max_estricto) {
descendencia <- forzar_n_max(
individuo = descendencia,
n_max = n_max_predictores
)
}
# Almacenar el nuevo descendiente en la nueva población: puede ocurrir que
# el individuo resultante del cruce y posterior mutación no contenga ningún
# predictor (todas sus posiciones son FALSE). Si esto ocurre, y para evitar
# que se produzca un error al ajustar el modelo, se introduce aleatoriamente
# un valor TRUE.
if (all(descendencia == FALSE)) {
descendencia[sample(x = 1:length(descendencia), size = 1)] <- TRUE
}
nueva_poblacion[j, ] <- descendencia
}
poblacion <- nueva_poblacion
# CRITERIO DE PARADA
# ==========================================================================
# Si durante las últimas n generaciones el mejor individuo no ha mejorado por
# una cantidad superior a tolerancia_parada, se detiene el algoritmo y no se
# crean nuevas generaciones.
if (parada_temprana && (i > rondas_parada)) {
ultimos_n <- tail(unlist(porcentaje_mejora), n = rondas_parada)
if (all(ultimos_n < tolerancia_parada)) {
print(paste(
"Algoritmo detenido en la generacion", i,
"por falta de mejora mínima de", tolerancia_parada,
"durante", rondas_parada,
"generaciones consecutivas."
))
break()
}
}
}
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------------
indice_mejor_individuo_global <- which.max(unlist(resultados_fitness))
mejor_fitness_global <- resultados_fitness[[indice_mejor_individuo_global]]
mejor_individuo_global <- resultados_individuo[[indice_mejor_individuo_global]]
# RESULTADOS
# ----------------------------------------------------------------------------
# La función devuelve una lista con 4 elementos:
# fitness:
# un vector con el fitness del mejor individuo de cada generación.
# mejor_individuo:
# una lista con la combinación de predictores del mejor individuo de cada
# generación.
# diferencia_abs:
# un vector con el porcentaje de mejora entre el mejor individuo de cada
# generación.
# porcentaje_mejora:
# un vector con el porcentaje de mejora entre el mejor individuo de cada
# generación.
# df_resultados:
# un dataframe con todos los resultados anteriores.
# mejor_individuo_global:
# la combinación de predictores del mejor individuo encontrado en todo el
# proceso.
resultados_fitness <- unlist(resultados_fitness)
diferencia_abs <- c(NA, unlist(diferencia_abs))
porcentaje_mejora <- c(NA, unlist(porcentaje_mejora))
# Para poder añadir al dataframe la secuencia de predictores, se concatenan.
predictores <- resultados_individuo[!sapply(resultados_individuo, is.null)]
predictores <- sapply(
X = predictores,
FUN = function(x) {
paste(x, collapse = ", ")
}
)
df_resultados <- data.frame(
generacion = seq_along(resultados_fitness),
fitness = resultados_fitness,
predictores = predictores,
porcentaje_mejora = porcentaje_mejora,
mejora = diferencia_abs
)
# Frecuencia relativa con la que aparece cada predictor en el mejor individuo
# a lo largo de las generaciones.
frecuencia_seleccion <- tibble::enframe(
x = table(unlist(resultados_individuo)),
name = "pred",
value = "frec"
)
frecuencia_seleccion <- frecuencia_seleccion[order(-frecuencia_seleccion$frec), ]
frecuencia_seleccion$frec <- as.numeric(100*(frecuencia_seleccion$frec/i))
top_5_predictores <- head(frecuencia_seleccion, 5)$pred
resultados <- list(
fitness = resultados_fitness,
mejor_individuo = resultados_individuo,
diferencia_abs = diferencia_abs,
porcentaje_mejora = porcentaje_mejora,
df_resultados = df_resultados,
mejor_individuo_global = mejor_individuo_global,
frecuencia_seleccion = frecuencia_seleccion,
top_5_predictores = top_5_predictores
)
end_time <- Sys.time()
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(resultados, "class") <- "seleccion_ga"
attr(resultados, 'fecha_creacion') <- end_time
attr(resultados, 'duracion_seleccion') <- difftime(end_time, start_time)
attr(resultados, 'generaciones') <- i
attr(resultados, 'predictores_seleccionados') <- mejor_individuo_global
attr(resultados, 'top_5_predictores') <- top_5_predictores
attr(resultados, 'fitness_optimo') <- mejor_fitness_global
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'modelo') <- modelo
attr(resultados, 'metrica') <- metrica
attr(resultados, 'elitismo') <- elitismo
attr(resultados, 'prob_mut') <- prob_mut
attr(resultados, 'metodo_seleccion') <- metodo_seleccion
attr(resultados, 'metodo_cruce') <- metodo_cruce
attr(resultados, 'parada_temprana') <- parada_temprana
attr(resultados, 'rondas_parada') <- rondas_parada
attr(resultados, 'tolerancia_parada') <- tolerancia_parada
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
cat("--------------------", "\n")
cat("Selección finalizada", "\n")
cat("--------------------", "\n")
cat("Fecha finalización:", as.character(Sys.time()), "\n")
cat("Duración selección: ")
print(difftime(end_time, start_time))
cat("Número de generaciones:", i, "\n")
cat("Predictores seleccionados:", mejor_individuo_global,"\n")
cat("Fitness final:", mejor_fitness_global, "\n")
cat("Top 5 predictores:", top_5_predictores,"\n")
cat("Modelo empleado:", modelo, "\n")
cat("Métrica empleada:", metrica, "\n")
cat("\n")
return(resultados)
}Comparación
Se compara el tiempo necesario para ejecutar el ejemplo del Titanic empleando paralelización y sin ella.
library(titanic)
library(dplyr)
# Se eliminan para este ejemplo 3 variables.
datos <- titanic_train %>% select(-Name, -Ticket, -Cabin)
# Se eliminan observaciones incompletas.
datos <- na.omit(datos)
# Se convierten a factor las variables de tipo character.
datos <- datos %>% mutate_if(is.character, as.factor)
datos <- datos %>% mutate(Survived = as.factor(Survived))Sin paralelización
library(tictoc)
tic()
seleccion <- seleccionar_predictores_ga(
x = datos[, -2],
y = datos[, 2],
n_poblacion = 50,
n_generaciones = 10,
n_max_predictores = 8,
n_min_predictores = 1,
elitismo = 0.01,
prob_mut = 0.01,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
modelo = "randomforest",
metrica = "accuracy",
nivel_referencia = "1",
cv = 5,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = FALSE,
verbose = FALSE
)## --------------------
## Selección finalizada
## --------------------
## Fecha finalización: 2020-02-01 19:00:31
## Duración selección: Time difference of 32.89964 secs
## Número de generaciones: 10
## Predictores seleccionados: Pclass Sex Age SibSp Parch Fare
## Fitness final: 0.8304442
## Top 5 predictores: Age Fare Parch Pclass Sex
## Modelo empleado: randomforest
## Métrica empleada: accuracytoc()## 33.045 sec elapsedCon paralelización
tic()
seleccion <- seleccionar_predictores_ga(
x = datos[, -2],
y = datos[, 2],
n_poblacion = 50,
n_generaciones = 10,
n_max_predictores = 8,
n_min_predictores = 1,
elitismo = 0.01,
prob_mut = 0.01,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
modelo = "randomforest",
metrica = "accuracy",
nivel_referencia = "1",
cv = 5,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = TRUE,
verbose = FALSE
)## --------------------
## Selección finalizada
## --------------------
## Fecha finalización: 2020-02-01 19:00:56
## Duración selección: Time difference of 24.85694 secs
## Número de generaciones: 10
## Predictores seleccionados: Pclass Sex Age SibSp Fare
## Fitness final: 0.833271
## Top 5 predictores: Age Fare Pclass Sex SibSp
## Modelo empleado: randomforest
## Métrica empleada: accuracytoc()## 24.861 sec elapsedPuede observarse que, para las mismas generaciones, la versión paralelizada tarda menos tiempo. Con ranger::ranger() la reducción de tiempo por paralelizar es menor que con randomForest::randomForest() ya que el primero ya paraleliza por defecto el ajuste de los modelos.
Implementación en R6
La programación de un algoritmo genético es un problema que se ajusta muy bien al paradigma de programación orientada a objetos (objetos que contienen tanto atributos como métodos). Este tipo de programación puede hacerse en R mediante la clase R6. En el siguiente link puede encontrarse una implementación del código anterior pero empleando clases R6.
Bibliografía
John McCall, Genetic algorithms for modelling and optimisation, Journal of Computational and Applied Mathematics, Volume 184, Issue 1, 2005
https://www.neuraldesigner.com/blog/genetic_algorithms_for_feature_selection
https://en.wikipedia.org/wiki/Crossover_(genetic_algorithm)
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.