graf_caso(data, contrib, 102, pred_base)
-
[1] "Predicción base: 0.18, Predicción: 0.85"[1] "Predicción base: 0.17, Predicción: 0.85"
-
-
@@ -873,23 +873,23 @@
Ejemplo: dia calibracion_gpos |> gt() |> fmt_number(where(is_double), decimals = 3)
-
-
@@ -1446,23 +1446,23 @@
preds_tbl |> head() |> gt() |> fmt_number(where(is_double), decimals = 3)
-
-
diff --git a/10-clasificacion-calibracion_files/figure-html/unnamed-chunk-4-1.png b/10-clasificacion-calibracion_files/figure-html/unnamed-chunk-4-1.png
index 6fbb77d3..483a45e0 100644
Binary files a/10-clasificacion-calibracion_files/figure-html/unnamed-chunk-4-1.png and b/10-clasificacion-calibracion_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/12-arboles_files/figure-html/unnamed-chunk-17-1.png b/12-arboles_files/figure-html/unnamed-chunk-17-1.png
index b3b512f7..58364f2d 100644
Binary files a/12-arboles_files/figure-html/unnamed-chunk-17-1.png and b/12-arboles_files/figure-html/unnamed-chunk-17-1.png differ
diff --git a/13-arboles-boosting.html b/13-arboles-boosting.html
index eb5d820a..0b798a65 100644
--- a/13-arboles-boosting.html
+++ b/13-arboles-boosting.html
@@ -935,7 +935,7 @@
system.time(modelo <- xgb.train(d_entrena, verbose = 1, nrounds = 10000, params = params))
system.time(modelo <- xgb.train(d_entrena, verbose = 1, nrounds = 10000, params = params)) user system elapsed
- 11.252 0.351 5.954 xgb.save(model = modelo, fname = "./cache/casas_boost.xgb")
diff --git a/14-interpretacion.html b/14-interpretacion.html
index 50cac15b..c3aa93fd 100644
--- a/14-interpretacion.html
+++ b/14-interpretacion.html
@@ -972,7 +972,7 @@
Ejemplo
} graf_caso(data, contrib, 10, pred_base)
-
[1] "Predicción base: 0.18, Predicción: 0.82"[1] "Predicción base: 0.17, Predicción: 0.87"
Ejemplo
La implementación oficial de R es lm, que en general tiene buen desempeño para datos que caben en memoria:
@@ -797,7 +797,7 @@Ejemplo
[[1]]
[,1]
-[1,] -0.22045916
+[1,] -0.22045918
[2,] 0.23149671
[3,] 0.09673478
@@ -817,11 +817,11 @@ Ejemplo
[[1]]
[,1]
[1,] 0.0079817399
-[2,] 0.0019486139
+[2,] 0.0019486138
[3,] 0.0005532746
[[2]]
-[1] 0.0003491883
+[1] 0.0003491882coef(lm.fit(cbind(1, x_ent), y_ent))
diff --git a/81-apendice-descenso_files/figure-html/unnamed-chunk-19-1.png b/81-apendice-descenso_files/figure-html/unnamed-chunk-19-1.png
index c79547b0..b2c97897 100644
Binary files a/81-apendice-descenso_files/figure-html/unnamed-chunk-19-1.png and b/81-apendice-descenso_files/figure-html/unnamed-chunk-19-1.png differ
diff --git a/82-apendice-descenso-estocastico.html b/82-apendice-descenso-estocastico.html
index 5cb55a96..93b4c852 100644
--- a/82-apendice-descenso-estocastico.html
+++ b/82-apendice-descenso-estocastico.html
@@ -798,13 +798,13 @@
get_weights(modelo)
get_weights(modelo)[[1]]
- [,1]
-[1,] -2.515035629
-[2,] -0.010236964
-[3,] -0.007203732
+ [,1]
+[1,] -2.64986229
+[2,] -0.05415167
+[3,] 0.04267785
[[2]]
-[1] -0.3029693Y verificamos que concuerda con la salida de lm:
diff --git a/82-apendice-descenso-estocastico_files/figure-html/unnamed-chunk-24-1.png b/82-apendice-descenso-estocastico_files/figure-html/unnamed-chunk-24-1.png index 16c4072e..f3509cd0 100644 Binary files a/82-apendice-descenso-estocastico_files/figure-html/unnamed-chunk-24-1.png and b/82-apendice-descenso-estocastico_files/figure-html/unnamed-chunk-24-1.png differ diff --git a/search.json b/search.json index 4d08c30a..941f1cfd 100644 --- a/search.json +++ b/search.json @@ -228,7 +228,7 @@ "href": "06-redes-neuronales-1.html#interacciones-en-redes-neuronales", "title": "6 Redes neuronales (intro)", "section": "6.2 Interacciones en redes neuronales", - "text": "6.2 Interacciones en redes neuronales\nEs posible capturar interacciones con redes neuronales. Consideremos el siguiente ejemplo simple:\n\nf <- function(x1, x2){\n 2 + 0.1* x1 + 0.1 * x2 + 10 * (x1 - 0.5) * (x2 - 0.5)\n}\ndat <- expand.grid(x1 = seq(0, 1, 0.05), x2 = seq(0, 1, 0.05))\ndat <- dat |> mutate(f = f(x1, x2))\nggplot(dat, aes(x=x1, y=x2)) + geom_tile(aes(fill=f))\n\n\n\n\nEsta función puede entenderse como un o exclusivo: la respuesta es alta sólo cuando \\(x_1\\) y \\(x_2\\) tienen valores opuestos (\\(x_1\\) grande pero \\(x_2\\) chica y viceversa).\nNo es posible modelar correctamente esta función mediante el modelo lineal (sin interacciones). Pero podemos incluir la interacción en el modelo lineal o intentar usar una red neuronal. Primero simulamos unos datos y probamos el modelo logístico con y sin interacciones:\n\nset.seed(322)\nn <- 2000\ndat_ent <- tibble(x1 = rbeta(n, 1, 1), x2 = rbeta(n, 1, 1)) |>\n mutate(f = f(x1, x2)) |>\n mutate(y = f + rnorm(n, 0, 0.1))\nmod_1 <- lm(y ~ x1 + x2, data = dat_ent)\nmod_1\n\n\nCall:\nlm(formula = y ~ x1 + x2, data = dat_ent)\n\nCoefficients:\n(Intercept) x1 x2 \n 1.8936 0.2046 0.2097 \n\n\nEl resultado del modelo lineal no es bueno:\n\ntibble(y_hat = fitted(mod_1), y = dat_ent$y) |> \n ggplot(aes(x = y_hat, y = y)) + geom_point(color = \"red\") +\n geom_abline() +\n coord_obs_pred()\n\n\n\n\nSin embargo, agregando una interacción lo mejoramos considerablemente (examina la raíz del error cuadrático medio, por ejemplo):\n\nmod_2 <- lm(y ~ x1 + x2 + x1:x2, data = dat_ent)\nmod_2\n\n\nCall:\nlm(formula = y ~ x1 + x2 + x1:x2, data = dat_ent)\n\nCoefficients:\n(Intercept) x1 x2 x1:x2 \n 4.499 -4.895 -4.885 9.964 \n\ntibble(y_hat = fitted(mod_2), y = dat_ent$y) |> \n ggplot(aes(x = y_hat, y = y)) + geom_point(color = \"red\") +\n geom_abline() +\n coord_obs_pred()\n\n\n\n\nObservese la gran diferencia de error entre los dos modelos (en este caso, el sobreajuste no es un problema).\nAhora consideramos qué red neuronal puede ser apropiada.\n\ntensorflow::set_random_seed(421) \nmod_inter <- keras_model_sequential()\nmod_inter |> \n layer_dense(units = 4, activation = \"sigmoid\",\n name = \"capa_intermedia\", input_shape = c(2)) |>\n layer_dense(units = 1, name = \"capa_final\",\n activation = \"linear\") \n\n\nmod_inter |> compile(loss = \"mse\", \n optimizer = optimizer_sgd(learning_rate = 0.3, momentum = 0.5))\nhistoria <- mod_inter |> \n fit(dat_ent |> select(x1, x2) |> as.matrix(), dat_ent$y,\n batch_size = 20,\n epochs = 100, verbose = 0)\n\nVerificamos que esta red captura la interacción:\n\npreds <- predict(mod_inter,\n dat |> select(x1, x2) |> as.matrix())\ndat <- dat |> mutate(f_red = preds)\nggplot(dat, aes(x = x1, y = x2)) + \n geom_tile(aes(fill = f_red))\n\n\n\n\n\npreds_ent <- predict(mod_inter, dat_ent |> select(x1, x2) |> as.matrix())\ntibble(pred = preds_ent[,1], f = dat_ent$y) |> \n ggplot(aes(x = pred, y = f)) +\n geom_point() +\n geom_abline(colour = \"red\") +\n coord_obs_pred()\n\n\n\n\nAunque podemos extraer los cálculos de la red ajustada, vamos a hacer el cálculo de la red a mano. La función feed forward es:\n\nbeta <- get_weights(mod_inter)\nfeed_fow <- function(beta, x){\n a <- h(t(beta[[1]]) %*% x + as.matrix(beta[[2]], 2, 1)) \n f <- t(beta[[3]]) %*% a + as.matrix(beta[[4]], 1, 1)\n f\n}\n\nObservación: ¿cómo funciona esta red? Consideremos la capa intermedia (3 unidades) para cuatro casos: \\((0,0), (0,1), (1,0), (1,1)\\):\n\nmat_entrada <- tibble(x_1 = c(0,0,1,1), x_2 = c(0,1,0,1)) |> as.matrix()\ncapa_1 <- keras_model(inputs = mod_inter$input,\n outputs = get_layer(mod_inter, \"capa_intermedia\")$output)\npred_mat <- predict(capa_1, mat_entrada) |> round(2)\nrownames(pred_mat) <- c(\"apagadas\", \"segunda\", \"primera\", \"ambas\")\npred_mat\n\n [,1] [,2] [,3] [,4]\napagadas 0.15 0.00 0.64 0.12\nsegunda 0.59 0.05 0.09 0.01\nprimera 0.01 0.05 0.07 0.60\nambas 0.05 0.56 0.00 0.06\n\n\nLos pesos de la última capa son:\n\nbeta[3:4]\n\n[[1]]\n [,1]\n[1,] -5.422091\n[2,] 4.781280\n[3,] 5.289693\n[4,] -5.083941\n\n[[2]]\n[1] 2.370601\n\n\nEjercicio: interpreta la red en términos de qué unidades están encendidas (valor cercano a 1) o apagadas (valor cercano a 0). ¿Puedes ajustar este modelo con dos tres unidades intermedias? Haz varias pruebas: ¿qué dificultades encuentras?" + "text": "6.2 Interacciones en redes neuronales\nEs posible capturar interacciones con redes neuronales. Consideremos el siguiente ejemplo simple:\n\nf <- function(x1, x2){\n 2 + 0.1* x1 + 0.1 * x2 + 10 * (x1 - 0.5) * (x2 - 0.5)\n}\ndat <- expand.grid(x1 = seq(0, 1, 0.05), x2 = seq(0, 1, 0.05))\ndat <- dat |> mutate(f = f(x1, x2))\nggplot(dat, aes(x=x1, y=x2)) + geom_tile(aes(fill=f))\n\n\n\n\nEsta función puede entenderse como un o exclusivo: la respuesta es alta sólo cuando \\(x_1\\) y \\(x_2\\) tienen valores opuestos (\\(x_1\\) grande pero \\(x_2\\) chica y viceversa).\nNo es posible modelar correctamente esta función mediante el modelo lineal (sin interacciones). Pero podemos incluir la interacción en el modelo lineal o intentar usar una red neuronal. Primero simulamos unos datos y probamos el modelo logístico con y sin interacciones:\n\nset.seed(322)\nn <- 2000\ndat_ent <- tibble(x1 = rbeta(n, 1, 1), x2 = rbeta(n, 1, 1)) |>\n mutate(f = f(x1, x2)) |>\n mutate(y = f + rnorm(n, 0, 0.1))\nmod_1 <- lm(y ~ x1 + x2, data = dat_ent)\nmod_1\n\n\nCall:\nlm(formula = y ~ x1 + x2, data = dat_ent)\n\nCoefficients:\n(Intercept) x1 x2 \n 1.8936 0.2046 0.2097 \n\n\nEl resultado del modelo lineal no es bueno:\n\ntibble(y_hat = fitted(mod_1), y = dat_ent$y) |> \n ggplot(aes(x = y_hat, y = y)) + geom_point(color = \"red\") +\n geom_abline() +\n coord_obs_pred()\n\n\n\n\nSin embargo, agregando una interacción lo mejoramos considerablemente (examina la raíz del error cuadrático medio, por ejemplo):\n\nmod_2 <- lm(y ~ x1 + x2 + x1:x2, data = dat_ent)\nmod_2\n\n\nCall:\nlm(formula = y ~ x1 + x2 + x1:x2, data = dat_ent)\n\nCoefficients:\n(Intercept) x1 x2 x1:x2 \n 4.499 -4.895 -4.885 9.964 \n\ntibble(y_hat = fitted(mod_2), y = dat_ent$y) |> \n ggplot(aes(x = y_hat, y = y)) + geom_point(color = \"red\") +\n geom_abline() +\n coord_obs_pred()\n\n\n\n\nObservese la gran diferencia de error entre los dos modelos (en este caso, el sobreajuste no es un problema).\nAhora consideramos qué red neuronal puede ser apropiada.\n\ntensorflow::set_random_seed(421) \nmod_inter <- keras_model_sequential()\nmod_inter |> \n layer_dense(units = 4, activation = \"sigmoid\",\n name = \"capa_intermedia\", input_shape = c(2)) |>\n layer_dense(units = 1, name = \"capa_final\",\n activation = \"linear\") \n\n\nmod_inter |> compile(loss = \"mse\", \n optimizer = optimizer_sgd(learning_rate = 0.3, momentum = 0.5))\nhistoria <- mod_inter |> \n fit(dat_ent |> select(x1, x2) |> as.matrix(), dat_ent$y,\n batch_size = 20,\n epochs = 100, verbose = 0)\n\nVerificamos que esta red captura la interacción:\n\npreds <- predict(mod_inter,\n dat |> select(x1, x2) |> as.matrix())\ndat <- dat |> mutate(f_red = preds)\nggplot(dat, aes(x = x1, y = x2)) + \n geom_tile(aes(fill = f_red))\n\n\n\n\n\npreds_ent <- predict(mod_inter, dat_ent |> select(x1, x2) |> as.matrix())\ntibble(pred = preds_ent[,1], f = dat_ent$y) |> \n ggplot(aes(x = pred, y = f)) +\n geom_point() +\n geom_abline(colour = \"red\") +\n coord_obs_pred()\n\n\n\n\nAunque podemos extraer los cálculos de la red ajustada, vamos a hacer el cálculo de la red a mano. La función feed forward es:\n\nbeta <- get_weights(mod_inter)\nfeed_fow <- function(beta, x){\n a <- h(t(beta[[1]]) %*% x + as.matrix(beta[[2]], 2, 1)) \n f <- t(beta[[3]]) %*% a + as.matrix(beta[[4]], 1, 1)\n f\n}\n\nObservación: ¿cómo funciona esta red? Consideremos la capa intermedia (3 unidades) para cuatro casos: \\((0,0), (0,1), (1,0), (1,1)\\):\n\nmat_entrada <- tibble(x_1 = c(0,0,1,1), x_2 = c(0,1,0,1)) |> as.matrix()\ncapa_1 <- keras_model(inputs = mod_inter$input,\n outputs = get_layer(mod_inter, \"capa_intermedia\")$output)\npred_mat <- predict(capa_1, mat_entrada) |> round(2)\nrownames(pred_mat) <- c(\"apagadas\", \"segunda\", \"primera\", \"ambas\")\npred_mat\n\n [,1] [,2] [,3] [,4]\napagadas 0.15 0.00 0.64 0.12\nsegunda 0.59 0.05 0.09 0.01\nprimera 0.01 0.05 0.07 0.60\nambas 0.05 0.56 0.00 0.06\n\n\nLos pesos de la última capa son:\n\nbeta[3:4]\n\n[[1]]\n [,1]\n[1,] -5.422077\n[2,] 4.781271\n[3,] 5.289690\n[4,] -5.083926\n\n[[2]]\n[1] 2.370607\n\n\nEjercicio: interpreta la red en términos de qué unidades están encendidas (valor cercano a 1) o apagadas (valor cercano a 0). ¿Puedes ajustar este modelo con dos tres unidades intermedias? Haz varias pruebas: ¿qué dificultades encuentras?" }, { "objectID": "06-redes-neuronales-1.html#cálculo-en-redes-feed-forward", @@ -263,7 +263,7 @@ "href": "06-redes-neuronales-1.html#ajuste-de-parámetros-introducción", "title": "6 Redes neuronales (intro)", "section": "6.7 Ajuste de parámetros (introducción)", - "text": "6.7 Ajuste de parámetros (introducción)\nConsideramos la versión con regularización ridge (también llamada L2) de la devianza de entrenamiento como nuestro función objetivo:\n\n\n\n\n\n\nAjuste de redes neuronales\n\n\n\nPara un problema de regresión, ajustamos los pesos \\(\\Theta^{(1)},\\Theta^{(2)},\\ldots, \\Theta^{(L)}\\) de la red intentando minimizar el error cuadrático medio (penalizado) sobre la muestra de entrenamiento: \\[D = -\\frac{1}{2n}\\sum_{i=1}^n (y^{(i)} - f(x^{i}))^2= + \\lambda \\sum_{l=2}^{L} \\sum_{k=1}^{n_{l-1}} \\sum_{j=1}^{n_l}(\\theta_{j,k}^{(l)})^2.\\] Este problema en general no es convexo y puede tener múltiples mínimos.\n\n\nVeremos el proceso de ajuste, selección de arquitectura, etc. más adelante. Por el momento hacemos unas observaciones acerca de este problema de minimización:\n\nHay varios algoritmos para minimizar este error, algunos avanzados incluyendo información de segundo orden (como Newton), pero actualmente las técnicas más populares, para redes grandes, están derivadas de descenso en gradiente. Más específicamente, una variación, que es descenso estocástico.\nQue el algoritmo depende principalmente de multiplicaciones de matrices y acumulaciones implica que puede escalarse de diversas maneras. Una es paralelizando sobre la muestra de entrenamiento (y calcular acumulados al final), pero también se pueden paralelizar las multiplicaciones de matrices (para lo cual los GPUs se prestan muy bien).\nPara redes neuronales, el gradiente se calcula con un algoritmo que se llama back-propagation, que es una aplicación de la regla de la cadena para propagar errores desde la capa de salida a lo largo de todas las capas para ajustar los pesos y sesgos.\nEn estos problemas no buscamos el mínimo global, sino un mínimo local de buen desempeño. Puede haber múltiples mínimos, puntos silla, regiones relativamente planas, precipicios (curvatura alta). Nótese que la simetría implica que podemos obtener la misma red cambiando pesos entre neuronas y las conexiones correspondientes. Esto implica que necesariamente hay varios mínimos.\nTodo esto dificulta el entrenamiento de redes neuronales grandes. Para redes grandes, ni siquiera esperamos a alcanzar un mínimo local, sino que nos a veces detenemos prematuramente cuando obtenemos el mejor desempeño posible.\nPara este problema, no tiene sentido comenzar las iteraciones con todos los pesos igual a cero, pues las unidades de la red son simétricas: no hay nada que distinga una de otra si todos los pesos son iguales. Esto quiere decir que si iteramos, todas las neuronas van a aprender lo mismo.\nEs importante no comenzar valores de los pesos grandes, pues las funciones logísticas pueden quedar en regiones planas donde la minimización es lenta, o podemos tener gradientes demasiado grandes y produzcan inestabilidad en el cálculo del gradiente.\nEl ajuste de la tasa de aprendizaje es un parámetro importante, más delicado que para problemas convexos. Generalmente lo tratamos con un hiperparámetro más que hay que afinar. Tasas demasiado grandes pueden llevarnos a mínimos locales relativamente malos.\nGeneralmente los pesos se inicializan al azar con variables independientes gaussianas o uniformes centradas en cero, y con varianza chica (por ejemplo \\(U(-0.5,0.5)\\)). Una recomendación es usar \\(U(-1/\\sqrt{m}, 1/\\sqrt{m})\\) donde \\(m\\) es el número de entradas. En general, hay que experimentar con este parámetro.\n\nEl proceso para ajustar una red es entonces:\n\nDefinir número de capas ocultas, número de neuronas por cada capa, y un valor del parámetro de regularización. Estandarizar las entradas. Usualmente podemos probar comenzar con una o dos capas ocultas, de tamaño proporcional al número de entradas. Es buena idea comenzar con una red relativamente grande que tienen error bajo de entrenamiento aunque sobreajuste, y después regularizar y refinar su tamaño.\nSeleccionar parámetros al azar para \\(\\Theta^{(2)},\\Theta^{(3)},\\ldots, \\Theta^{(L)}\\). Se toman, por ejemplo, normales con media 0 y varianza chica.\nCorrer un algoritmo de minimización del error mostrada arriba. Es necesario experimentar con los parámetros del algoritmo de minimización.\nVerificar convergencia del algoritmo a un mínimo local (o el algoritmo no está mejorando).\nPredecir usando el modelo ajustado.\n\nFinalmente, podemos probar distintas arquitecturas y valores del parámetros de regularización, para afinar estos parámetros según validación cruzada o una muestra de validación.\n\nEjemplo (regresión)\n\ndat_grasa <- read_csv(file = '../datos/bodyfat.csv') \nset.seed(183)\ngrasa_particion <- initial_split(dat_grasa, 0.5)\ngrasa_ent <- training(grasa_particion)\ngrasa_pr <- testing(grasa_particion)\nnrow(grasa_ent)\n\n[1] 126\n\n\nUna exploración de este conjunto de datos revela algunos datos sospechosos. En particular un individuo con estatura de 30 pulgadas (alrededor de 75 cm), con peso normal. Probablemente no queremos incluir en entrenamiento este caso, y tampoco hacer predicciones para posibles personas que tengan tales dimensiones:\n\nlibrary(patchwork)\ng_1 <- grasa_ent |> ggplot(aes(x = estatura, y = peso)) + \n geom_point()\ng_2 <- grasa_ent |> ggplot(aes(x = abdomen, y = peso)) + \n geom_point()\ng_1 + g_2\n\n\n\n\n\ngrasa_receta <- recipe(grasacorp ~ ., grasa_ent) |> \n step_filter(estatura > 50) |>\n step_normalize(all_predictors()) |> \n prep()\n\n\nlibrary(keras)\n# entrenamiento\nx_grasa <- grasa_receta |> juice() |> \n select(-grasacorp) |> as.matrix()\nvars_nombres <- colnames(x_grasa)\ny_grasa <- grasa_receta |> juice() |> pull(grasacorp)\n# validación\nx_grasa_pr <- grasa_receta |> bake(grasa_pr) |> \n select(-grasacorp) |> as.matrix()\ny_grasa_pr <- grasa_receta |> bake(grasa_pr) |> pull(grasacorp)\n\n\nmodelo_red <- keras_model_sequential() |> \n layer_dense(units = 50, activation = \"sigmoid\") |> \n layer_dense(units = 50, activation = \"sigmoid\") |> \n layer_dense(units = 1, activation = \"linear\")\nmodelo_red |> compile(\n loss = \"mse\", metrics = \"mse\",\n optimizer = optimizer_sgd(learning_rate = 0.01, momentum = 0.9)\n)\n# esto es más eficiente hacerlo con callbacks en general:\nhistoria <- modelo_red |> fit(\n x = x_grasa, y = y_grasa,\n validation_data = list(x_grasa_pr, y_grasa_pr),\n batch_size = 30, epochs = 250, verbose = 1)\n\nVemos que con esta red podemos alcanzar un error de entrenamiento cercano a cero, aún cuando vemos que sobreajusta al evaluar con la muestra de validación.\n\nplot(historia, smooth = FALSE)\n\n\n\n\nNotamos que las predicciones no son muy buenas:\n\npreds <- predict(modelo_red, x_grasa_pr) \npreds |> head()\n\n [,1]\n[1,] 9.117618\n[2,] 22.663441\n[3,] 8.569652\n[4,] 8.558598\n[5,] 8.842549\n[6,] 8.559493\n\n\nY obtenemos el siguiente resultado:\n\ng_1 <- tibble(preds = preds[, 1], y = y_grasa_pr) |> \n ggplot(aes(x = preds, y = y)) + \n geom_point() + \n geom_abline(slope = 1, intercept = 0, color = \"red\") +\n coord_obs_pred()\ng_1\n\n\n\n\nPodemos ahora experimentar con los parámetros del optimizador, número de unidades, número de capas y regularización L2.\n\nCada cambio de número de unidades/capas o regularización requiere ajustes a la tasa de aprendizaje y otros parámetros del optimizador.\nVarias arquitecturas (número de capas y unidades) pueden dar resultados similares. En este caso, usualmente escogemos el modelo más computacionalmente simple, o dependiendo del tipo de errores de cada modelo.\nRecordamos que una consecuencia del sobreajuste es que el error de prueba es mayor que el de entrenamiento (hay un gap de generalización). Para mejorar el desempeño, podemos reducir el error de entrenamiento (reducir sesgo), y/o reducir sobreajuste (gap entre entrenamiento y prueba).\n\n\nmodelo_red_2 <- keras_model_sequential() |> \n layer_dense(units = 30, activation = \"sigmoid\", \n kernel_regularizer = regularizer_l2(0.1)) |> \n layer_dense(units = 30, activation = \"sigmoid\",\n kernel_regularizer = regularizer_l2(0.1)) |> \n layer_dense(units = 1, activation = \"linear\", \n kernel_regularizer = regularizer_l2(0.01))\nmodelo_red_2 |> compile(loss = \"mse\",metrics = c(\"mse\"),\n optimizer = optimizer_sgd(learning_rate = 0.0005, momentum = 0.95)\n)\n# esto es más eficiente hacerlo con callbacks en general:\nhistoria <- modelo_red_2 |> fit(\n x = x_grasa, y = y_grasa,\n validation_data = list(x_grasa_pr, y_grasa_pr),\n batch_size = 30, epochs = 500, verbose = 1)\n\n\nplot(historia, smooth = FALSE)\n\n\n\n\nLa afinación nos da mejores resultados:\n\npreds_2 <- predict(modelo_red_2, x_grasa_pr) \ng_2 <- tibble(preds_2 = preds_2[, 1], y = y_grasa_pr) |> \n ggplot(aes(x = preds_2, y = y)) + \n geom_point() + \n geom_abline(slope = 1, intercept = 0, color = \"red\") +\n coord_obs_pred()\ng_1 + g_2\n\n\n\n\n\n\nResumen\nEjercicio: En nuestra referencia The Deep Learning Book, puedes revisar la sección 11.4.1 (afinación manual) y 11.4.2 (afinación por búsqueda automática). Explica qué efecto tienen sobre sesgo y sobreajuste los siguientes parámetros:\n\nNúmero de capas\nNúmero de unidades por capa oculta\nRegularización L2 de los pesos\nTasa de aprendizaje (ojo, este parámetro es especial)\n\nEn la parte 2 de redes neuronales veremos más estrategias de regularización, y cómo agregar estructura a redes adaptada a problemas particulares. Notamos para por el momento que las redes neuronales con solamente capas conexas no han sido particularmente exitosas para ningún problema particular. Típicamente requieren una gran cantidad de datos y entrenarlas es más difícil que otros métodos. Sin embargo, la historia es diferente cuando escogemos adecuadamente la arquitectura de la red para un problema dado, es decir, qué unidades están conectadas y bajo qué patrones. En estos casos, las conexiones pueden ser relativamente pocas comparadas con redes totalmente conexas con el mismo tamaño de unidades, y si la estructura es correcta, no incurrimos en mucho sesgo. Ejemplos son procesamiento de imágenes, audio, o lenguaje natural." + "text": "6.7 Ajuste de parámetros (introducción)\nConsideramos la versión con regularización ridge (también llamada L2) de la devianza de entrenamiento como nuestro función objetivo:\n\n\n\n\n\n\nAjuste de redes neuronales\n\n\n\nPara un problema de regresión, ajustamos los pesos \\(\\Theta^{(1)},\\Theta^{(2)},\\ldots, \\Theta^{(L)}\\) de la red intentando minimizar el error cuadrático medio (penalizado) sobre la muestra de entrenamiento: \\[D = -\\frac{1}{2n}\\sum_{i=1}^n (y^{(i)} - f(x^{i}))^2= + \\lambda \\sum_{l=2}^{L} \\sum_{k=1}^{n_{l-1}} \\sum_{j=1}^{n_l}(\\theta_{j,k}^{(l)})^2.\\] Este problema en general no es convexo y puede tener múltiples mínimos.\n\n\nVeremos el proceso de ajuste, selección de arquitectura, etc. más adelante. Por el momento hacemos unas observaciones acerca de este problema de minimización:\n\nHay varios algoritmos para minimizar este error, algunos avanzados incluyendo información de segundo orden (como Newton), pero actualmente las técnicas más populares, para redes grandes, están derivadas de descenso en gradiente. Más específicamente, una variación, que es descenso estocástico.\nQue el algoritmo depende principalmente de multiplicaciones de matrices y acumulaciones implica que puede escalarse de diversas maneras. Una es paralelizando sobre la muestra de entrenamiento (y calcular acumulados al final), pero también se pueden paralelizar las multiplicaciones de matrices (para lo cual los GPUs se prestan muy bien).\nPara redes neuronales, el gradiente se calcula con un algoritmo que se llama back-propagation, que es una aplicación de la regla de la cadena para propagar errores desde la capa de salida a lo largo de todas las capas para ajustar los pesos y sesgos.\nEn estos problemas no buscamos el mínimo global, sino un mínimo local de buen desempeño. Puede haber múltiples mínimos, puntos silla, regiones relativamente planas, precipicios (curvatura alta). Nótese que la simetría implica que podemos obtener la misma red cambiando pesos entre neuronas y las conexiones correspondientes. Esto implica que necesariamente hay varios mínimos.\nTodo esto dificulta el entrenamiento de redes neuronales grandes. Para redes grandes, ni siquiera esperamos a alcanzar un mínimo local, sino que nos a veces detenemos prematuramente cuando obtenemos el mejor desempeño posible.\nPara este problema, no tiene sentido comenzar las iteraciones con todos los pesos igual a cero, pues las unidades de la red son simétricas: no hay nada que distinga una de otra si todos los pesos son iguales. Esto quiere decir que si iteramos, todas las neuronas van a aprender lo mismo.\nEs importante no comenzar valores de los pesos grandes, pues las funciones logísticas pueden quedar en regiones planas donde la minimización es lenta, o podemos tener gradientes demasiado grandes y produzcan inestabilidad en el cálculo del gradiente.\nEl ajuste de la tasa de aprendizaje es un parámetro importante, más delicado que para problemas convexos. Generalmente lo tratamos con un hiperparámetro más que hay que afinar. Tasas demasiado grandes pueden llevarnos a mínimos locales relativamente malos.\nGeneralmente los pesos se inicializan al azar con variables independientes gaussianas o uniformes centradas en cero, y con varianza chica (por ejemplo \\(U(-0.5,0.5)\\)). Una recomendación es usar \\(U(-1/\\sqrt{m}, 1/\\sqrt{m})\\) donde \\(m\\) es el número de entradas. En general, hay que experimentar con este parámetro.\n\nEl proceso para ajustar una red es entonces:\n\nDefinir número de capas ocultas, número de neuronas por cada capa, y un valor del parámetro de regularización. Estandarizar las entradas. Usualmente podemos probar comenzar con una o dos capas ocultas, de tamaño proporcional al número de entradas. Es buena idea comenzar con una red relativamente grande que tienen error bajo de entrenamiento aunque sobreajuste, y después regularizar y refinar su tamaño.\nSeleccionar parámetros al azar para \\(\\Theta^{(2)},\\Theta^{(3)},\\ldots, \\Theta^{(L)}\\). Se toman, por ejemplo, normales con media 0 y varianza chica.\nCorrer un algoritmo de minimización del error mostrada arriba. Es necesario experimentar con los parámetros del algoritmo de minimización.\nVerificar convergencia del algoritmo a un mínimo local (o el algoritmo no está mejorando).\nPredecir usando el modelo ajustado.\n\nFinalmente, podemos probar distintas arquitecturas y valores del parámetros de regularización, para afinar estos parámetros según validación cruzada o una muestra de validación.\n\nEjemplo (regresión)\n\ndat_grasa <- read_csv(file = '../datos/bodyfat.csv') \nset.seed(183)\ngrasa_particion <- initial_split(dat_grasa, 0.5)\ngrasa_ent <- training(grasa_particion)\ngrasa_pr <- testing(grasa_particion)\nnrow(grasa_ent)\n\n[1] 126\n\n\nUna exploración de este conjunto de datos revela algunos datos sospechosos. En particular un individuo con estatura de 30 pulgadas (alrededor de 75 cm), con peso normal. Probablemente no queremos incluir en entrenamiento este caso, y tampoco hacer predicciones para posibles personas que tengan tales dimensiones:\n\nlibrary(patchwork)\ng_1 <- grasa_ent |> ggplot(aes(x = estatura, y = peso)) + \n geom_point()\ng_2 <- grasa_ent |> ggplot(aes(x = abdomen, y = peso)) + \n geom_point()\ng_1 + g_2\n\n\n\n\n\ngrasa_receta <- recipe(grasacorp ~ ., grasa_ent) |> \n step_filter(estatura > 50) |>\n step_normalize(all_predictors()) |> \n prep()\n\n\nlibrary(keras)\n# entrenamiento\nx_grasa <- grasa_receta |> juice() |> \n select(-grasacorp) |> as.matrix()\nvars_nombres <- colnames(x_grasa)\ny_grasa <- grasa_receta |> juice() |> pull(grasacorp)\n# validación\nx_grasa_pr <- grasa_receta |> bake(grasa_pr) |> \n select(-grasacorp) |> as.matrix()\ny_grasa_pr <- grasa_receta |> bake(grasa_pr) |> pull(grasacorp)\n\n\nmodelo_red <- keras_model_sequential() |> \n layer_dense(units = 50, activation = \"sigmoid\") |> \n layer_dense(units = 50, activation = \"sigmoid\") |> \n layer_dense(units = 1, activation = \"linear\")\nmodelo_red |> compile(\n loss = \"mse\", metrics = \"mse\",\n optimizer = optimizer_sgd(learning_rate = 0.01, momentum = 0.9)\n)\n# esto es más eficiente hacerlo con callbacks en general:\nhistoria <- modelo_red |> fit(\n x = x_grasa, y = y_grasa,\n validation_data = list(x_grasa_pr, y_grasa_pr),\n batch_size = 30, epochs = 250, verbose = 1)\n\nVemos que con esta red podemos alcanzar un error de entrenamiento cercano a cero, aún cuando vemos que sobreajusta al evaluar con la muestra de validación.\n\nplot(historia, smooth = FALSE)\n\n\n\n\nNotamos que las predicciones no son muy buenas:\n\npreds <- predict(modelo_red, x_grasa_pr) \npreds |> head()\n\n [,1]\n[1,] 9.117623\n[2,] 22.663464\n[3,] 8.569628\n[4,] 8.558584\n[5,] 8.842519\n[6,] 8.559466\n\n\nY obtenemos el siguiente resultado:\n\ng_1 <- tibble(preds = preds[, 1], y = y_grasa_pr) |> \n ggplot(aes(x = preds, y = y)) + \n geom_point() + \n geom_abline(slope = 1, intercept = 0, color = \"red\") +\n coord_obs_pred()\ng_1\n\n\n\n\nPodemos ahora experimentar con los parámetros del optimizador, número de unidades, número de capas y regularización L2.\n\nCada cambio de número de unidades/capas o regularización requiere ajustes a la tasa de aprendizaje y otros parámetros del optimizador.\nVarias arquitecturas (número de capas y unidades) pueden dar resultados similares. En este caso, usualmente escogemos el modelo más computacionalmente simple, o dependiendo del tipo de errores de cada modelo.\nRecordamos que una consecuencia del sobreajuste es que el error de prueba es mayor que el de entrenamiento (hay un gap de generalización). Para mejorar el desempeño, podemos reducir el error de entrenamiento (reducir sesgo), y/o reducir sobreajuste (gap entre entrenamiento y prueba).\n\n\nmodelo_red_2 <- keras_model_sequential() |> \n layer_dense(units = 30, activation = \"sigmoid\", \n kernel_regularizer = regularizer_l2(0.1)) |> \n layer_dense(units = 30, activation = \"sigmoid\",\n kernel_regularizer = regularizer_l2(0.1)) |> \n layer_dense(units = 1, activation = \"linear\", \n kernel_regularizer = regularizer_l2(0.01))\nmodelo_red_2 |> compile(loss = \"mse\",metrics = c(\"mse\"),\n optimizer = optimizer_sgd(learning_rate = 0.0005, momentum = 0.95)\n)\n# esto es más eficiente hacerlo con callbacks en general:\nhistoria <- modelo_red_2 |> fit(\n x = x_grasa, y = y_grasa,\n validation_data = list(x_grasa_pr, y_grasa_pr),\n batch_size = 30, epochs = 500, verbose = 1)\n\n\nplot(historia, smooth = FALSE)\n\n\n\n\nLa afinación nos da mejores resultados:\n\npreds_2 <- predict(modelo_red_2, x_grasa_pr) \ng_2 <- tibble(preds_2 = preds_2[, 1], y = y_grasa_pr) |> \n ggplot(aes(x = preds_2, y = y)) + \n geom_point() + \n geom_abline(slope = 1, intercept = 0, color = \"red\") +\n coord_obs_pred()\ng_1 + g_2\n\n\n\n\n\n\nResumen\nEjercicio: En nuestra referencia The Deep Learning Book, puedes revisar la sección 11.4.1 (afinación manual) y 11.4.2 (afinación por búsqueda automática). Explica qué efecto tienen sobre sesgo y sobreajuste los siguientes parámetros:\n\nNúmero de capas\nNúmero de unidades por capa oculta\nRegularización L2 de los pesos\nTasa de aprendizaje (ojo, este parámetro es especial)\n\nEn la parte 2 de redes neuronales veremos más estrategias de regularización, y cómo agregar estructura a redes adaptada a problemas particulares. Notamos para por el momento que las redes neuronales con solamente capas conexas no han sido particularmente exitosas para ningún problema particular. Típicamente requieren una gran cantidad de datos y entrenarlas es más difícil que otros métodos. Sin embargo, la historia es diferente cuando escogemos adecuadamente la arquitectura de la red para un problema dado, es decir, qué unidades están conectadas y bajo qué patrones. En estos casos, las conexiones pueden ser relativamente pocas comparadas con redes totalmente conexas con el mismo tamaño de unidades, y si la estructura es correcta, no incurrimos en mucho sesgo. Ejemplos son procesamiento de imágenes, audio, o lenguaje natural." }, { "objectID": "07-intervalos-predictivos.html#inferencia-predictiva-conforme", @@ -361,7 +361,7 @@ "href": "08-clasificacion-1.html#ejercicio-datos-de-diabetes", "title": "8 Clasificación y probabilidad", "section": "8.8 Ejercicio: datos de diabetes", - "text": "8.8 Ejercicio: datos de diabetes\nYa están divididos los datos en entrenamiento y prueba\n\ndiabetes_ent <- as_tibble(MASS::Pima.tr)\ndiabetes_pr <- as_tibble(MASS::Pima.te)\ndiabetes_ent |> head() |> gt()\n\n\n\n\n\n \n \n \n npreg\n glu\n bp\n skin\n bmi\n ped\n age\n type\n \n \n \n 5\n86\n68\n28\n30.2\n0.364\n24\nNo\n 7\n195\n70\n33\n25.1\n0.163\n55\nYes\n 5\n77\n82\n41\n35.8\n0.156\n35\nNo\n 0\n165\n76\n43\n47.9\n0.259\n26\nNo\n 0\n107\n60\n25\n26.4\n0.133\n23\nNo\n 5\n97\n76\n27\n35.6\n0.378\n52\nYes\n \n \n \n\n\n\ndiabetes_ent$id <- 1:nrow(diabetes_ent)\ndiabetes_pr$id <- 1:nrow(diabetes_pr)\n\nAunque no es necesario, podemos normalizar:\n\nreceta_diabetes <- recipe(type ~ ., diabetes_ent) |>\n update_role(id, new_role = \"id_variable\") |> \n step_normalize(all_numeric()) \ndiabetes_ent_s <- receta_diabetes |> prep() |> juice() \ndiabetes_pr_s <- receta_diabetes |> prep() |> bake(diabetes_pr)\n\n\nmodelo_lineal <- logistic_reg(mode = \"classification\") |> \n set_engine(\"glm\")\nflujo_diabetes <- workflow() |> \n add_model(modelo_lineal) |> \n add_recipe(receta_diabetes)\nflujo_ajustado <- fit(flujo_diabetes, diabetes_ent)\nsaveRDS(flujo_ajustado, \"cache/flujo_ajustado_diabetes.rds\")\nflujo_ajustado\n\n══ Workflow [trained] ══════════════════════════════════════════════════════════\nPreprocessor: Recipe\nModel: logistic_reg()\n\n── Preprocessor ────────────────────────────────────────────────────────────────\n1 Recipe Step\n\n• step_normalize()\n\n── Model ───────────────────────────────────────────────────────────────────────\n\nCall: stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)\n\nCoefficients:\n(Intercept) npreg glu bp skin bmi \n -0.95583 0.34734 1.01705 -0.05473 -0.02247 0.51263 \n ped age \n 0.55928 0.45201 \n\nDegrees of Freedom: 199 Total (i.e. Null); 192 Residual\nNull Deviance: 256.4 \nResidual Deviance: 178.4 AIC: 194.4\n\n\nAhora calculamos devianza de prueba y error de clasificación:\n\npreds_prueba <- \n predict(flujo_ajustado, diabetes_pr, type= \"prob\") |> \n bind_cols(predict(flujo_ajustado, diabetes_pr)) |> \n bind_cols(diabetes_pr |> select(type))\npreds_prueba\n\n# A tibble: 332 × 4\n .pred_No .pred_Yes .pred_class type \n <dbl> <dbl> <fct> <fct>\n 1 0.232 0.768 Yes Yes \n 2 0.960 0.0403 No No \n 3 0.975 0.0253 No No \n 4 0.959 0.0413 No Yes \n 5 0.204 0.796 Yes Yes \n 6 0.265 0.735 Yes Yes \n 7 0.590 0.410 No Yes \n 8 0.780 0.220 No No \n 9 0.558 0.442 No No \n10 0.798 0.202 No Yes \n# ℹ 322 more rows\n\n\n\nlevels(preds_prueba$type)\n\n[1] \"No\" \"Yes\"\n\n# ponemos event_level si \"positivo\" no es el primer factor\nmetricas <- metric_set(accuracy, mn_log_loss)\nmetricas(preds_prueba, truth = type, .pred_Yes, estimate = .pred_class, \n event_level = \"second\")\n\n# A tibble: 2 × 3\n .metric .estimator .estimate\n <chr> <chr> <dbl>\n1 accuracy binary 0.801\n2 mn_log_loss binary 0.441\n\n\nVamos a repetir usando keras.\n\nlibrary(keras)\nx_ent <- diabetes_ent_s |> select(-type, -id) |> as.matrix()\ny_ent <- diabetes_ent_s$type == \"Yes\"\nx_prueba <- diabetes_pr_s |> select(-type, -id) |> as.matrix()\ny_prueba <- diabetes_pr_s$type == 'Yes'\n# definición de estructura del modelo (regresión logística)\n# es posible hacerlo con workflows como vimos arriba, \n# pero aquí usamos directamente la interfaz de keras en R\nn_entrena <- nrow(x_ent)\nmodelo_diabetes <- keras_model_sequential() |>\n layer_dense(units = 1, #una sola respuesta,\n activation = \"sigmoid\", # combinar variables linealmente y aplicar función logística\n kernel_initializer = initializer_constant(0), #inicializamos coeficientes en 0\n bias_initializer = initializer_constant(0)) #inicializamos ordenada en 0\n# compilar seleccionando cantidad a minimizar, optimizador y métricas\nmodelo_diabetes |> compile(\n loss = \"binary_crossentropy\", # devianza es entropía cruzada\n optimizer = optimizer_sgd(learning_rate = 0.75), # descenso en gradiente\n metrics = list(\"binary_crossentropy\"))\n# Ahora iteramos\n# Primero probamos con un número bajo de iteraciones\nhistoria <- modelo_diabetes |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 10 # número de iteraciones\n)\nplot(historia)\n\n\n\n\nY ahora podemos correr más iteraciones adicionales:\n\nhistoria <- modelo_diabetes |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 1000, # número de iteraciones\n verbose = 0\n)\n\nLos errores de entrenamiento y prueba son:\n\nevaluate(modelo_diabetes, x_ent, y_ent)\n\n loss binary_crossentropy \n 0.4459766 0.4459766 \n\n\n\nevaluate(modelo_diabetes, x_prueba, y_prueba)\n\n loss binary_crossentropy \n 0.4406986 0.4406986 \n\n\nVeamos que coeficientes obtuvimos:\n\nget_weights(modelo_diabetes)\n\n[[1]]\n [,1]\n[1,] 0.34734303\n[2,] 1.01705003\n[3,] -0.05472932\n[4,] -0.02247143\n[5,] 0.51263177\n[6,] 0.55927491\n[7,] 0.45200703\n\n[[2]]\n[1] -0.9558301\n\n\nque coinciden con los valores que obtuvimos usando regresión logística de glm. La única diferencia es que el algoritmo de optimización que se usa en cada caso es diferente: con keras utilizamos descenso en gradiente, mientras que glm usa Newton-Raphson.\n\nflujo_ajustado |> extract_fit_parsnip()\n\nparsnip model object\n\n\nCall: stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)\n\nCoefficients:\n(Intercept) npreg glu bp skin bmi \n -0.95583 0.34734 1.01705 -0.05473 -0.02247 0.51263 \n ped age \n 0.55928 0.45201 \n\nDegrees of Freedom: 199 Total (i.e. Null); 192 Residual\nNull Deviance: 256.4 \nResidual Deviance: 178.4 AIC: 194.4" + "text": "8.8 Ejercicio: datos de diabetes\nYa están divididos los datos en entrenamiento y prueba\n\ndiabetes_ent <- as_tibble(MASS::Pima.tr)\ndiabetes_pr <- as_tibble(MASS::Pima.te)\ndiabetes_ent |> head() |> gt()\n\n\n\n\n\n \n \n \n npreg\n glu\n bp\n skin\n bmi\n ped\n age\n type\n \n \n \n 5\n86\n68\n28\n30.2\n0.364\n24\nNo\n 7\n195\n70\n33\n25.1\n0.163\n55\nYes\n 5\n77\n82\n41\n35.8\n0.156\n35\nNo\n 0\n165\n76\n43\n47.9\n0.259\n26\nNo\n 0\n107\n60\n25\n26.4\n0.133\n23\nNo\n 5\n97\n76\n27\n35.6\n0.378\n52\nYes\n \n \n \n\n\n\ndiabetes_ent$id <- 1:nrow(diabetes_ent)\ndiabetes_pr$id <- 1:nrow(diabetes_pr)\n\nAunque no es necesario, podemos normalizar:\n\nreceta_diabetes <- recipe(type ~ ., diabetes_ent) |>\n update_role(id, new_role = \"id_variable\") |> \n step_normalize(all_numeric()) \ndiabetes_ent_s <- receta_diabetes |> prep() |> juice() \ndiabetes_pr_s <- receta_diabetes |> prep() |> bake(diabetes_pr)\n\n\nmodelo_lineal <- logistic_reg(mode = \"classification\") |> \n set_engine(\"glm\")\nflujo_diabetes <- workflow() |> \n add_model(modelo_lineal) |> \n add_recipe(receta_diabetes)\nflujo_ajustado <- fit(flujo_diabetes, diabetes_ent)\nsaveRDS(flujo_ajustado, \"cache/flujo_ajustado_diabetes.rds\")\nflujo_ajustado\n\n══ Workflow [trained] ══════════════════════════════════════════════════════════\nPreprocessor: Recipe\nModel: logistic_reg()\n\n── Preprocessor ────────────────────────────────────────────────────────────────\n1 Recipe Step\n\n• step_normalize()\n\n── Model ───────────────────────────────────────────────────────────────────────\n\nCall: stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)\n\nCoefficients:\n(Intercept) npreg glu bp skin bmi \n -0.95583 0.34734 1.01705 -0.05473 -0.02247 0.51263 \n ped age \n 0.55928 0.45201 \n\nDegrees of Freedom: 199 Total (i.e. Null); 192 Residual\nNull Deviance: 256.4 \nResidual Deviance: 178.4 AIC: 194.4\n\n\nAhora calculamos devianza de prueba y error de clasificación:\n\npreds_prueba <- \n predict(flujo_ajustado, diabetes_pr, type= \"prob\") |> \n bind_cols(predict(flujo_ajustado, diabetes_pr)) |> \n bind_cols(diabetes_pr |> select(type))\npreds_prueba\n\n# A tibble: 332 × 4\n .pred_No .pred_Yes .pred_class type \n <dbl> <dbl> <fct> <fct>\n 1 0.232 0.768 Yes Yes \n 2 0.960 0.0403 No No \n 3 0.975 0.0253 No No \n 4 0.959 0.0413 No Yes \n 5 0.204 0.796 Yes Yes \n 6 0.265 0.735 Yes Yes \n 7 0.590 0.410 No Yes \n 8 0.780 0.220 No No \n 9 0.558 0.442 No No \n10 0.798 0.202 No Yes \n# ℹ 322 more rows\n\n\n\nlevels(preds_prueba$type)\n\n[1] \"No\" \"Yes\"\n\n# ponemos event_level si \"positivo\" no es el primer factor\nmetricas <- metric_set(accuracy, mn_log_loss)\nmetricas(preds_prueba, truth = type, .pred_Yes, estimate = .pred_class, \n event_level = \"second\")\n\n# A tibble: 2 × 3\n .metric .estimator .estimate\n <chr> <chr> <dbl>\n1 accuracy binary 0.801\n2 mn_log_loss binary 0.441\n\n\nVamos a repetir usando keras.\n\nlibrary(keras)\nx_ent <- diabetes_ent_s |> select(-type, -id) |> as.matrix()\ny_ent <- diabetes_ent_s$type == \"Yes\"\nx_prueba <- diabetes_pr_s |> select(-type, -id) |> as.matrix()\ny_prueba <- diabetes_pr_s$type == 'Yes'\n# definición de estructura del modelo (regresión logística)\n# es posible hacerlo con workflows como vimos arriba, \n# pero aquí usamos directamente la interfaz de keras en R\nn_entrena <- nrow(x_ent)\nmodelo_diabetes <- keras_model_sequential() |>\n layer_dense(units = 1, #una sola respuesta,\n activation = \"sigmoid\", # combinar variables linealmente y aplicar función logística\n kernel_initializer = initializer_constant(0), #inicializamos coeficientes en 0\n bias_initializer = initializer_constant(0)) #inicializamos ordenada en 0\n# compilar seleccionando cantidad a minimizar, optimizador y métricas\nmodelo_diabetes |> compile(\n loss = \"binary_crossentropy\", # devianza es entropía cruzada\n optimizer = optimizer_sgd(learning_rate = 0.75), # descenso en gradiente\n metrics = list(\"binary_crossentropy\"))\n# Ahora iteramos\n# Primero probamos con un número bajo de iteraciones\nhistoria <- modelo_diabetes |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 10 # número de iteraciones\n)\nplot(historia)\n\n\n\n\nY ahora podemos correr más iteraciones adicionales:\n\nhistoria <- modelo_diabetes |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 1000, # número de iteraciones\n verbose = 0\n)\n\nLos errores de entrenamiento y prueba son:\n\nevaluate(modelo_diabetes, x_ent, y_ent)\n\n loss binary_crossentropy \n 0.4459766 0.4459766 \n\n\n\nevaluate(modelo_diabetes, x_prueba, y_prueba)\n\n loss binary_crossentropy \n 0.4406986 0.4406986 \n\n\nVeamos que coeficientes obtuvimos:\n\nget_weights(modelo_diabetes)\n\n[[1]]\n [,1]\n[1,] 0.34734306\n[2,] 1.01705003\n[3,] -0.05472932\n[4,] -0.02247145\n[5,] 0.51263183\n[6,] 0.55927497\n[7,] 0.45200703\n\n[[2]]\n[1] -0.9558302\n\n\nque coinciden con los valores que obtuvimos usando regresión logística de glm. La única diferencia es que el algoritmo de optimización que se usa en cada caso es diferente: con keras utilizamos descenso en gradiente, mientras que glm usa Newton-Raphson.\n\nflujo_ajustado |> extract_fit_parsnip()\n\nparsnip model object\n\n\nCall: stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)\n\nCoefficients:\n(Intercept) npreg glu bp skin bmi \n -0.95583 0.34734 1.01705 -0.05473 -0.02247 0.51263 \n ped age \n 0.55928 0.45201 \n\nDegrees of Freedom: 199 Total (i.e. Null); 192 Residual\nNull Deviance: 256.4 \nResidual Deviance: 178.4 AIC: 194.4" }, { "objectID": "08-clasificacion-1.html#probabilidades-y-pérdida-0-1", @@ -417,7 +417,7 @@ "href": "09-clasificacion-2.html#perfil-de-un-clasificador-binario-y-curvas-roc", "title": "9 Decisiones de clasificación", "section": "9.5 Perfil de un clasificador binario y curvas ROC", - "text": "9.5 Perfil de un clasificador binario y curvas ROC\nEn lugar de examinar cada punto de corte por separado, podemos hacer el análisis de todos los posibles puntos de corte mediante la curva ROC (receiver operating characteristic, de ingeniería).\n\n\n\n\n\n\nTip\n\n\n\nPara un problema de clasificación binaria, dadas estimaciones \\(\\hat{p}_1(x)\\), la curva ROC grafica todos los pares de (1-especificidad, sensibilidad) para cada posible punto de corte \\(\\hat{p}_1(x) > \\alpha\\).\n\n\n\nEjemplo\n\nroc_tbl <- roc_curve(preds_prueba, \n truth = type, .pred_Yes)\nroc_tbl\n\n# A tibble: 321 × 3\n .threshold specificity sensitivity\n <dbl> <dbl> <dbl>\n 1 -Inf 0 1\n 2 0.00305 0 1\n 3 0.00345 0.00448 1\n 4 0.00406 0.00897 1\n 5 0.00437 0.0135 1\n 6 0.00475 0.0224 1\n 7 0.00485 0.0269 1\n 8 0.00504 0.0314 1\n 9 0.00512 0.0359 1\n10 0.00562 0.0404 1\n# ℹ 311 more rows\n\nggplot(roc_tbl, \n aes(x = 1 - specificity, y = sensitivity)) +\n geom_path(aes(colour = .threshold), size = 1.2) +\n geom_abline(colour = \"gray\") + \n coord_equal() +\n xlab(\"Tasa de falsos positivos\") + \n ylab(\"Sensibilidad\")\n\nWarning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.\nℹ Please use `linewidth` instead.\n\n\n\n\n\nEn esta gráfica podemos ver todos los clasificadores posibles basados en las probabilidades de clase. Podemos usar estas curvas como evaluación de nuestros clasificadores, dejando para más tarde la selección del punto de corte, si esto es necesario (por ejemplo, dependiendo de los costos de cada tipo de error).\nTambién podemos definir una medida resumen del desempeño de un clasificador según esta curva:\n\n\n\n\n\n\nTip\n\n\n\nLa medida AUC (area under the curve) para un clasificador es el área bajo la curva generada por los pares sensibilidad-especificidad de la curva ROC.\n\n\n\nroc_auc(preds_prueba, type, .pred_Yes)\n\n# A tibble: 1 × 3\n .metric .estimator .estimate\n <chr> <chr> <dbl>\n1 roc_auc binary 0.761\n\n\nCuanto más cerca de uno es este valor, mejor discriminan las probabilidades.\nTambién es útil para comparar modelos. Consideremos el modelo de los datos de diabetes que incluyen todas las variables:\n\nmod_2 <- logistic_reg() |> \n set_engine(\"keras\") |> \n set_mode(\"classification\") |>\n set_args(epochs = 250, \n optimizer = optimizer_sgd(lr = 0.3),\n batch_size = nrow(diabetes_ent), \n verbose = FALSE) |> \n fit(type ~ ., juice(receta_diabetes))\npreds_prueba_completo <- predict(\n mod_2, prueba_baked, type ='prob') |> \n bind_cols(prueba_baked)\npreds_prueba_2 <- bind_rows(\n preds_prueba |> mutate(modelo = \"IMC y edad\"),\n preds_prueba_completo |> mutate(modelo = \"Todas\")) |> \n group_by(modelo)\n\nY graficamos juntas:\n\nroc_2_tbl <- roc_curve(preds_prueba_2, type, .pred_Yes)\nggplot(roc_2_tbl, aes(x = 1 - specificity, y = sensitivity)) +\n geom_path(aes(colour = modelo), size = 1.2) +\n geom_abline(colour = \"gray\") + \n coord_equal() +\n xlab(\"Tasa de falsos positivos\") + \n ylab(\"Sensibilidad\")\n\n\n\n\nComparación auc:\n\nroc_auc(preds_prueba_2, type, .pred_Yes) \n\n# A tibble: 2 × 4\n modelo .metric .estimator .estimate\n <chr> <chr> <chr> <dbl>\n1 IMC y edad roc_auc binary 0.761\n2 Todas roc_auc binary 0.866\n\n\nEn este ejemplo, vemos que casi no importa qué perfil de especificidad y sensibilidad busquemos: el clasificador que usa todas las variables domina casi siempre al clasificador que sólo utiliza IMC y edad." + "text": "9.5 Perfil de un clasificador binario y curvas ROC\nEn lugar de examinar cada punto de corte por separado, podemos hacer el análisis de todos los posibles puntos de corte mediante la curva ROC (receiver operating characteristic, de ingeniería).\n\n\n\n\n\n\nTip\n\n\n\nPara un problema de clasificación binaria, dadas estimaciones \\(\\hat{p}_1(x)\\), la curva ROC grafica todos los pares de (1-especificidad, sensibilidad) para cada posible punto de corte \\(\\hat{p}_1(x) > \\alpha\\).\n\n\n\nEjemplo\n\nroc_tbl <- roc_curve(preds_prueba, \n truth = type, .pred_Yes)\nroc_tbl\n\n# A tibble: 321 × 3\n .threshold specificity sensitivity\n <dbl> <dbl> <dbl>\n 1 -Inf 0 1\n 2 0.00305 0 1\n 3 0.00345 0.00448 1\n 4 0.00406 0.00897 1\n 5 0.00437 0.0135 1\n 6 0.00476 0.0224 1\n 7 0.00486 0.0269 1\n 8 0.00504 0.0314 1\n 9 0.00512 0.0359 1\n10 0.00562 0.0404 1\n# ℹ 311 more rows\n\nggplot(roc_tbl, \n aes(x = 1 - specificity, y = sensitivity)) +\n geom_path(aes(colour = .threshold), size = 1.2) +\n geom_abline(colour = \"gray\") + \n coord_equal() +\n xlab(\"Tasa de falsos positivos\") + \n ylab(\"Sensibilidad\")\n\nWarning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.\nℹ Please use `linewidth` instead.\n\n\n\n\n\nEn esta gráfica podemos ver todos los clasificadores posibles basados en las probabilidades de clase. Podemos usar estas curvas como evaluación de nuestros clasificadores, dejando para más tarde la selección del punto de corte, si esto es necesario (por ejemplo, dependiendo de los costos de cada tipo de error).\nTambién podemos definir una medida resumen del desempeño de un clasificador según esta curva:\n\n\n\n\n\n\nTip\n\n\n\nLa medida AUC (area under the curve) para un clasificador es el área bajo la curva generada por los pares sensibilidad-especificidad de la curva ROC.\n\n\n\nroc_auc(preds_prueba, type, .pred_Yes)\n\n# A tibble: 1 × 3\n .metric .estimator .estimate\n <chr> <chr> <dbl>\n1 roc_auc binary 0.761\n\n\nCuanto más cerca de uno es este valor, mejor discriminan las probabilidades.\nTambién es útil para comparar modelos. Consideremos el modelo de los datos de diabetes que incluyen todas las variables:\n\nmod_2 <- logistic_reg() |> \n set_engine(\"keras\") |> \n set_mode(\"classification\") |>\n set_args(epochs = 250, \n optimizer = optimizer_sgd(lr = 0.3),\n batch_size = nrow(diabetes_ent), \n verbose = FALSE) |> \n fit(type ~ ., juice(receta_diabetes))\npreds_prueba_completo <- predict(\n mod_2, prueba_baked, type ='prob') |> \n bind_cols(prueba_baked)\npreds_prueba_2 <- bind_rows(\n preds_prueba |> mutate(modelo = \"IMC y edad\"),\n preds_prueba_completo |> mutate(modelo = \"Todas\")) |> \n group_by(modelo)\n\nY graficamos juntas:\n\nroc_2_tbl <- roc_curve(preds_prueba_2, type, .pred_Yes)\nggplot(roc_2_tbl, aes(x = 1 - specificity, y = sensitivity)) +\n geom_path(aes(colour = modelo), size = 1.2) +\n geom_abline(colour = \"gray\") + \n coord_equal() +\n xlab(\"Tasa de falsos positivos\") + \n ylab(\"Sensibilidad\")\n\n\n\n\nComparación auc:\n\nroc_auc(preds_prueba_2, type, .pred_Yes) \n\n# A tibble: 2 × 4\n modelo .metric .estimator .estimate\n <chr> <chr> <chr> <dbl>\n1 IMC y edad roc_auc binary 0.761\n2 Todas roc_auc binary 0.864\n\n\nEn este ejemplo, vemos que casi no importa qué perfil de especificidad y sensibilidad busquemos: el clasificador que usa todas las variables domina casi siempre al clasificador que sólo utiliza IMC y edad." }, { "objectID": "09-clasificacion-2.html#curvas-de-precisión-sensibilidad", @@ -613,7 +613,7 @@ "href": "13-arboles-boosting.html#submuestreo", "title": "13 Métodos basados en árboles: boosting", "section": "13.11 Submuestreo", - "text": "13.11 Submuestreo\nPuede funcionar bien construir cada uno de los árboles con submuestras de la muestra de entrenamiento, como una manera adicional de reducir varianza al construir nuestro predictor (esta idea es parecida a la de los bosques aleatorios, aquí igualmente perturbamos la muestra de entrenamiento en cada paso para evitar sobreajuste). Adicionalmente, este proceso acelera considerablemente las iteraciones de boosting, y en algunos casos sin penalización en desempeño.\nEn boosting se pueden tomar submuestras (una fracción mayor a 0.5 de la muestra de entrenamiento, pero puede ser más chica para conjuntos grandes de entrenamiento) sin reemplazo.\nEste parámetro también puede ser afinado con muestra de validación o validación cruzada.\nFinalmente, hacemos una evaluación correcta de validación cruzada:\n\n#xgboost es el default\nmodelo_boosting <- boost_tree(learn_rate = 0.003, trees = tune(), \n mtry = 5, tree_depth = 4) |> \n set_mode(\"regression\") |> \n set_args(objective = \"reg:squarederror\")\nflujo_casas <- workflow() |> add_recipe(receta_casas) |> add_model(modelo_boosting)\n\n\nnum_arboles_tbl <- tibble(trees = seq(1, 10000, 100))\nset.seed(81)\nparticion_vc <- vfold_cv(casas_entrena, v = 10)\nmis_metricas <- metric_set(mape, rsq)\nresultados <- tune_grid(flujo_casas, particion_vc, grid = num_arboles_tbl, metrics = mis_metricas)\n\n\ncollect_metrics(resultados) |> \n filter(.metric == \"mape\", trees > 10) |> \n ggplot(aes(x = trees, y = mean, ymin = mean - std_err, ymax = mean + std_err)) + \n geom_point() +\n geom_line() + geom_ribbon(alpha = 0.2) + ylab(\"mape val_cruzada\")\n\n\n\n\n\ncollect_metrics(resultados) |> filter(.metric == \"mape\") |> arrange(mean) |> head()\n\n# A tibble: 6 × 7\n trees .metric .estimator mean n std_err .config \n <dbl> <chr> <chr> <dbl> <int> <dbl> <chr> \n1 9901 mape standard 9.62 10 0.516 Preprocessor1_Model100\n2 9801 mape standard 9.62 10 0.516 Preprocessor1_Model099\n3 9701 mape standard 9.62 10 0.516 Preprocessor1_Model098\n4 9301 mape standard 9.62 10 0.515 Preprocessor1_Model094\n5 9401 mape standard 9.62 10 0.515 Preprocessor1_Model095\n6 9201 mape standard 9.62 10 0.515 Preprocessor1_Model093\n\n\nFinalmente podemos guardar el modelo en un formato estándar (R, Python, GCP y otros):\n\nsystem.time(modelo <- xgb.train(d_entrena, verbose = 1, nrounds = 10000, params = params))\n\n user system elapsed \n 11.252 0.351 5.954 \n\nxgb.save(model = modelo, fname = \"./cache/casas_boost.xgb\")\n\n[1] TRUE" + "text": "13.11 Submuestreo\nPuede funcionar bien construir cada uno de los árboles con submuestras de la muestra de entrenamiento, como una manera adicional de reducir varianza al construir nuestro predictor (esta idea es parecida a la de los bosques aleatorios, aquí igualmente perturbamos la muestra de entrenamiento en cada paso para evitar sobreajuste). Adicionalmente, este proceso acelera considerablemente las iteraciones de boosting, y en algunos casos sin penalización en desempeño.\nEn boosting se pueden tomar submuestras (una fracción mayor a 0.5 de la muestra de entrenamiento, pero puede ser más chica para conjuntos grandes de entrenamiento) sin reemplazo.\nEste parámetro también puede ser afinado con muestra de validación o validación cruzada.\nFinalmente, hacemos una evaluación correcta de validación cruzada:\n\n#xgboost es el default\nmodelo_boosting <- boost_tree(learn_rate = 0.003, trees = tune(), \n mtry = 5, tree_depth = 4) |> \n set_mode(\"regression\") |> \n set_args(objective = \"reg:squarederror\")\nflujo_casas <- workflow() |> add_recipe(receta_casas) |> add_model(modelo_boosting)\n\n\nnum_arboles_tbl <- tibble(trees = seq(1, 10000, 100))\nset.seed(81)\nparticion_vc <- vfold_cv(casas_entrena, v = 10)\nmis_metricas <- metric_set(mape, rsq)\nresultados <- tune_grid(flujo_casas, particion_vc, grid = num_arboles_tbl, metrics = mis_metricas)\n\n\ncollect_metrics(resultados) |> \n filter(.metric == \"mape\", trees > 10) |> \n ggplot(aes(x = trees, y = mean, ymin = mean - std_err, ymax = mean + std_err)) + \n geom_point() +\n geom_line() + geom_ribbon(alpha = 0.2) + ylab(\"mape val_cruzada\")\n\n\n\n\n\ncollect_metrics(resultados) |> filter(.metric == \"mape\") |> arrange(mean) |> head()\n\n# A tibble: 6 × 7\n trees .metric .estimator mean n std_err .config \n <dbl> <chr> <chr> <dbl> <int> <dbl> <chr> \n1 9901 mape standard 9.62 10 0.516 Preprocessor1_Model100\n2 9801 mape standard 9.62 10 0.516 Preprocessor1_Model099\n3 9701 mape standard 9.62 10 0.516 Preprocessor1_Model098\n4 9301 mape standard 9.62 10 0.515 Preprocessor1_Model094\n5 9401 mape standard 9.62 10 0.515 Preprocessor1_Model095\n6 9201 mape standard 9.62 10 0.515 Preprocessor1_Model093\n\n\nFinalmente podemos guardar el modelo en un formato estándar (R, Python, GCP y otros):\n\nsystem.time(modelo <- xgb.train(d_entrena, verbose = 1, nrounds = 10000, params = params))\n\n user system elapsed \n 9.440 0.265 5.062 \n\nxgb.save(model = modelo, fname = \"./cache/casas_boost.xgb\")\n\n[1] TRUE" }, { "objectID": "13-arboles-boosting.html#otros-parámetros", @@ -669,7 +669,7 @@ "href": "14-interpretacion.html#explicación-de-predicciones", "title": "14 Interpretación de modelos", "section": "14.5 Explicación de predicciones", - "text": "14.5 Explicación de predicciones\nUna tarea común que se puede requerir para transparentar las predicciones es distribuir contribuciones a la predicción de cada variable.\nConsideremos una predictor \\(f(X)\\), y hacemos una predicción \\(f(x)\\) para un caso \\(x = (x_1,x_2,\\ldots, x_p)\\) particular. ¿Cómo contribuye cada variable a la predicción f(x)?\nUna manera común (ojo: no es la única) de asignar contribuciones si el modelo es lineal y no tiene interacciones es la siguiente:\nSi \\(f(X_1,\\ldots, X_p) = \\beta_0 + \\beta_1 X_1 + \\cdots + \\beta_p X_p\\) , y \\(x^* = (x_1^*,x_2^*,\\ldots, x_p^*)\\) es un caso o instancia particular, podemos definir la contribución \\(\\phi_j = \\phi_j (x^*)\\) de la variable \\(j\\) en la predicción como\n\\[\\phi_j = \\beta_j(x_j^* - \\mu_j),\\] donde \\(\\mu_j\\) es la media de la variable \\(x_j\\) (por ejemplo en el conjunto de entrenamiento). Podemos también definir fácilmente la contribución de un subconjunto \\(W\\) de variables dada, como \\[\\phi_W (x) = \\sum_{j\\in W} \\beta_j(x_j^* - \\mu_j)\\].\nContribución es un buen nombre para estas cantidades, pues satisfacen (usando la linealidad y las definiciones de arriba):\n\nLas contribuciones suman la diferencia de la predicción con la predicción media: \\[\\sum_j \\phi_j(x) = f(x) - E(f(x))\\]\nSi una variable satisface \\(\\phi_W (x) = \\phi_{W\\cup j} (x)\\) para cualquier subconjunto de variables \\(W\\), entonces \\(\\phi_j\\) es cero.\nSi nuestro predictor se puede escribir como \\(f(x) = f_1(x) + f_2(x)\\) entonces la contribución de cada variable en la predicción \\(f(x)\\) es la suma de las contribución en \\(f_1\\) y en \\(f_2\\)\nPara cualquier subconjunto de variables, podemos considerar la contribución de una variable \\(j\\) cuando agregamos al subcojunto la variable \\(j\\) como \\(\\phi_{W \\cup j} (x) - \\phi_W (x)\\). Si para cualquier subconjunto \\(W\\) tenemos que las variables \\(j\\) y \\(k\\) tienen la misma contribución al agregarlos a \\(W\\) (para toda \\(W\\)), entonces \\(\\phi_j =\\phi_k\\).\n\nAhora la pregunta es: ¿cómo definimos las contribuciones para un predictor \\(f\\) más complejo? (por ejemplo, un bosque aleatorio). Resulta que el concepto de contribución o atribución es uno que se estudia en teoría de juegos, y se puede demostrar que hay una sola forma de hacerlo cumpliendo las propiedades señaladas arriba: mediante los valores de Shapley.\nAntes de seguir, veremos las dificultades que podemos encontrar para definir la atribución o contribuciones:\n\nEjemplo\nConsideramos un modelo lineal para empezar. Por ejemplo \\[f(x_1, x_2) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2,\\] ¿Cómo podemos atribuir contribuciones de \\(x_1\\) y \\(x_2\\) para una entrada particular \\(x^* = (x_1^*, x_2^*).\\)\nEn primer lugar, podríamos ver qué pasa cuando pensamos que no tenemos ninguna variable disponible. Esta cantidad la podríamos definir como (promediando sobre la muestra de entrenamiento): \\[v(0) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{(i)}) = \\beta_0 + \\beta_1 \\bar{x_1} + \\beta_2\\bar{x_2}\\]\nque es la predicción media (usamos linealidad de \\(f\\)).\nAhora supongamos que usamos la variable \\(x_1\\). En este caso, podríamos calcular una predicción promediando sobre la variable \\(x_2\\), que no tenemos disponible:\n\\[v(1) = \\frac{1}{N}\\sum_i f(x_1^{*},x_2^{(i)}) = \\beta_0 + \\beta_1 x_1^* + \\beta_2\\bar{x_2}\\] Entonces la contribución cuando agregamos la variable 1 es: \\[v(1) -v(0) = \\beta_1(x_1^* - \\bar{x_1})\\] Ahora supongamos que ya tenemos la variable 2, y agregamos la variable 1. La predicción con todas las variables es \\[v(1,2) = \\beta_0 + \\beta_1 x_1^* + \\beta_2 x_2^*\\] y cuando solo tenemos la variable \\(2\\) es \\[v(2) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{*}) = \\beta_0 + \\beta_1\\bar{x_1} + \\beta_2 x_2^*, \\] y la contribución en este caso es: \\[v(1,2) - v(2) = \\beta_1 (x_1^* - \\bar{x_1}),\\] que es la misma cantidad que \\(v(1)-v(0)\\). Así que no importa como “incluyamos” la variable 1 bajo este criterio, el resultado es el mismo, y podemos definir la contribución de \\(x_1\\) como definimos arriba para un predictor lineal.\nAhora consideremos el caso más complicado donde tenemos una interacción multiplicativa:\n\\[f(x_1, x_2) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_3x_1x_2\\]\nVeremos qué pasa cuando seguimos la misma idea de arriba: en primer lugar, tenemos\n\\[v(0) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{(i)}) =\n\\beta_0 + \\beta_1 \\overline{x_1} + \\beta_2\\overline{x_2} + \\beta_3\\overline{x_1x_2}\\]\ny \\[v(1) = \\frac{1}{N}\\sum_i f(x_1^{*},x_2^{(i)}) = \\beta_0 + \\beta_1 x_1^* + \\beta_2\\bar{x_2} + (\\beta_3 \\bar{x_2})x_1^*\\] La contribución sería \\[v(1)-v(0) = \\beta_1(x_1^* - \\bar{x_1}) +\\beta_3\\bar{x_2}\\left(x_1^*-\\frac{\\overline{x_1x_2}}{\\bar{x_2}}\\right)\\] Por otro lado, igual que arriba\n\\[v(2) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{*}) = \\beta_0 + \\beta_1\\bar{x_1} + \\beta_2 x_2^* + \\beta_3x_2^*\\bar{x_1}\\]\ny entonces: \\[v(1,2) - v(2) = \\beta_1(x_1^*-\\bar{x_1})+ \\beta_3x_2^*(x_1^* - \\bar{x_1})\\] Y las cantidades \\(v(1) - v(0)\\) y \\(v(1,2)- v(2)\\) no son iguales en este caso. Tenemos dos órdenes distintos para poner las variables en el modelo, así que ponemos:\n\\[\\phi_1 = \\frac{1}{2}(v(1,2)-v(2)) + \\frac{1}{2}(v(1)-v(0)).\\] y análogamente para \\[\\phi_2 = \\frac{1}{2}(v(1,2)-v(1)) + \\frac{1}{2}(v(2)-v(0)).\\]\nNótese que estas cantidades satisfacen: \\(\\phi_1 + \\phi_2 = v(1,2) -v(0) = f(x_1^*, x_2^*) - \\frac{1}{N}\\sum_i f(x_1^{(i)}, x_2^{(i)}),\\)\nes decir, la diferencia entre la predicción que nos interesa y la predicción media, que es la primera propiedad que buscábamos.\nConsideremos ahora el caso de 3 variables, y nos interesa encontrar la contribución de la variable 1. Ahora tenemos 6 órdenes distintos para poner las variables en el predictor. Dos de ellos terminan con \\(v(1,2,3) - v(2,3)\\), uno de ellos contiene la comparación \\(v(1,2)-v(2)\\), uno de ellos contiene la comparación \\(v(1,3)-v(3)\\) y dos de ellos contiene la comparación \\(v(1)-v(0)\\). Así que ponderaríamos:\n\\[\\phi_1 = \\frac{2}{6}(v(1,2,3) - v(2,3)) + \\frac{1}{6}(v(1,2)-v(2)) + \\frac{1}{6}(v(1,3)-v(3)) + \\frac{2}{6}(v(1)-v(0))\\] Puedes checar que \\(\\phi_1+ \\phi_2 +\\phi_3 = f(x_1^*, x_2^*, x_3^*) - \\frac{1}{N}\\sum_i f(x_1^{(i)}, x_2^{(i)}, x_3^{(i)})\\). Los valores así definidos satisfacen todas las propiedades que discutimos arriba.\n\n\n14.5.1 Valores de Shapley\nLa idea principal en los valores de Shapley es que para entender la contribución de una variable, es necesario considerar su contribución según distintos órdenes en el que consideramos las variables.\nDefinimos primero la contribución de un conjunto de variables \\(S\\subset M,\\) donde \\(M =\\{x_1, x_2\\ldots, x_p\\}\\), marginalizando como sigue con la muestra de entrenamiento:\nSuponiendo que \\(S = \\{1,2,3,\\ldots ,k\\}\\) entonces \\[v(S) = \\frac{1}{N}\\sum_i f(x_1^*,x_2^*,\\ldots, x_k^*, x_k^{(i)}, x_{k+1}^{(i)}, \\ldots, x_p^{(i)})\\] Donde se mantienen fijos los elementos de \\(S\\), y promediamos sobre los elementos que no están en \\(S\\). Esto podemos hacerlo para cualquier subconjunto \\(S\\) de variables. Ahora consideramos una variable dada \\(j\\), y consideramos las diferencias: \\[v(S\\cup {j}) - v(S)\\] Es decir, cuánto contribuye la variable \\(j\\) cuando la agregamos a las variables que ya teníamos en \\(S\\).\nEl valor de Shapley de la variable \\(j\\) es\n\\[\\phi_j = \\sum_{S\\subset M, j\\notin S} \\frac{|S|!(M - |S| -1)!}{M!} (v(S\\cup j) -v(S))\\]\nLos ponderadores están dados como discutimos arriba: hay \\(M!\\) posibles ordenamientos de las variables. Los ordenamientos que incluyen la comparación \\(v(S\\cup j) -v(S)\\) satisfacen: las primeras \\(|S|\\) variables pueden haber entrado de \\(|S|!\\) maneras, luego sigue la \\(j\\), y el restante de las variables pueden entrar de \\((M-|S|-1)!\\) maneras.\n\n\nEjemplo: precios de casas\nConsideremos entonces el bosque que construimos arriba para los precios de casas, y consideramos una casa particular:\n\nset.seed(232)\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[21,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\n\nLa suma de los valores phi da (aproximadamente, por el algoritmo usado) la diferencia entre nuestra predicción para esta instancia y la predicción promedio:\n\nsum(shapley_val$results$phi)\n\n[1] -10.40783\n\n\nAhora hacemos otro ejemplo:\n\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[52,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\nsum(shapley_val$results$phi)\n\n[1] -47.69727\n\n\n\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[121,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\n\nObservaciones:\n\nNótese que las contribuciones de cada variable en un mismo nivel puede ser diferente en diferentes instancias, pues los modelos que usamos típicamente incluyen interacciones.\nCalcular sobre todas las posibles permutaciones de las variables es demasiado costoso. Para estimar el valor de Shapley, en iml se toma una muestra de permutaciones, y para cada una se calcula el valor correspondiente \\(v(S\\cup j) - v(S),\\) dependiendo dónde aparece \\(j\\), y promediando sobre las variables que no aparecen en \\(S\\) como mostramos arriba. No es necesario ponderar la muestra de permutaciones.\nExisten también versiones adaptadas a los árboles que son más rápìdas.\n\n\n\nEjemplo\nPara un modelo de clasificación, xgboost calcula los valores de shapley en la escala logit (del predictor lineal)\n\nlibrary(xgboost)\nnse_tbl <- nse_entrena %>% ungroup %>% \n dplyr::select(ab, conex_inte, ocupados, \n num_tvd, educa_jefe, tam_loc, combustible, num_auto)\nmat_x <- model.matrix(ab ~ -1 + ., nse_tbl)\nd_entrena <- xgb.DMatrix(data = mat_x, label = nse_tbl$ab)\nboost_nse <- xgboost(mat_x, label = as.numeric(nse_tbl$ab == \"1\"), \n nrounds = 500, eta = 0.1, \n max_depth = 3, subsample = .5,\n objective = \"binary:logistic\", nthread = 2, verbose = 0)\nshap_nse <- xgb.plot.shap(mat_x, model = boost_nse, top_n = 8, n_col = 2)\n\n\n\npred_base <- mean(nse_tbl$ab==1)\n\n\ndata <- shap_nse$data\ncontrib <- shap_nse$shap_contrib\ndat_g <- tibble(educa_jefe = data[, \"educa_jefe\"], shap_educa = contrib[, \"educa_jefe\"],\n shap_ocupados = contrib[, \"ocupados\"], ) %>% \n sample_n(2000)\nggplot(dat_g, aes(x = educa_jefe, y = shap_educa, colour = shap_ocupados)) + \n geom_jitter(alpha = 0.9, width = 0.5, height = 0) +\n scale_color_gradient(low=\"orange\", high=\"black\")\n\n\n\n\nY podemos ver predicciones individuales:\n\ngraf_caso <- function(data, contrib, ind_caso, pred_base){\n dat_tbl <- tibble(variable = colnames(data), valor = data[ind_caso, ], \n shap = contrib[ind_caso,]) %>% \n unite(\"variable_valor\", variable, valor) %>% \n mutate(variable_valor = fct_reorder(variable_valor, shap))\n pred_logit <- log(pred_base / (1-pred_base)) + sum(dat_tbl$shap)\n pred <- 1 / (1 + exp(-pred_logit))\n sprintf(\"Predicción base: %0.2f, Predicción: %0.2f\", pred_base, pred) %>% print\n ggplot(dat_tbl, aes(x= variable_valor, y = shap, ymax = shap, ymin = 0)) + \n coord_flip() +geom_point() + geom_hline(yintercept = 0, colour = \"gray\") +\n geom_linerange() \n}\ngraf_caso(data, contrib, 10, pred_base)\n\n[1] \"Predicción base: 0.18, Predicción: 0.82\"\n\n\n\n\n\n\ngraf_caso(data, contrib, 102, pred_base)\n\n[1] \"Predicción base: 0.18, Predicción: 0.85\"\n\n\n\n\n\nDiscusión:\n\nUn valor de contribución que puede ser más apropiado para los valores Shapley es condicionando en cada caso a la información que se tiene durante el proceso de adición de variables. Es decir, usamos \\[v(S) = E_{X_C|X_S}\\left [ f(X_S,X_C) \\, | \\, X_S = x_s^* \\right].\\] Esta cantidad es teóricamente más apropiada para hacer predicciones cuando “no tenemos” las variables de \\(X_C\\). Sin embargo, calcular esta cantidad es considerablemente difícil, pues requiere modelar también la conjunta de \\((X_1,\\ldots, X_p)\\), lo cuál en general es difícil. Hasta en regresión lineal, incluso sin interacciones, no es trivial hacer estos cálculos. Generalmente se utilizan como mostramos arriba. Una desventaja clara del proceso como lo mostramos arriba es que puede ser que al hacer estos promedios, usemos partes del modelo con poca información y malas predicciones. Los valores de Shapley pueden ser ruidosos en este caso.\nLos que mostramos arriba son llamados comunmente valores SHAP o explicación aditiva de Shapley. Existen muchas variaciones de la aplicación de valores de Shapley para entender predictores y es un tema de investigación activo.\n\n\n\n\n\nHastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2017. The Elements of Statistical Learning. Springer Series en Statistics. Springer New York Inc. http://web.stanford.edu/~hastie/ElemStatLearn/." + "text": "14.5 Explicación de predicciones\nUna tarea común que se puede requerir para transparentar las predicciones es distribuir contribuciones a la predicción de cada variable.\nConsideremos una predictor \\(f(X)\\), y hacemos una predicción \\(f(x)\\) para un caso \\(x = (x_1,x_2,\\ldots, x_p)\\) particular. ¿Cómo contribuye cada variable a la predicción f(x)?\nUna manera común (ojo: no es la única) de asignar contribuciones si el modelo es lineal y no tiene interacciones es la siguiente:\nSi \\(f(X_1,\\ldots, X_p) = \\beta_0 + \\beta_1 X_1 + \\cdots + \\beta_p X_p\\) , y \\(x^* = (x_1^*,x_2^*,\\ldots, x_p^*)\\) es un caso o instancia particular, podemos definir la contribución \\(\\phi_j = \\phi_j (x^*)\\) de la variable \\(j\\) en la predicción como\n\\[\\phi_j = \\beta_j(x_j^* - \\mu_j),\\] donde \\(\\mu_j\\) es la media de la variable \\(x_j\\) (por ejemplo en el conjunto de entrenamiento). Podemos también definir fácilmente la contribución de un subconjunto \\(W\\) de variables dada, como \\[\\phi_W (x) = \\sum_{j\\in W} \\beta_j(x_j^* - \\mu_j)\\].\nContribución es un buen nombre para estas cantidades, pues satisfacen (usando la linealidad y las definiciones de arriba):\n\nLas contribuciones suman la diferencia de la predicción con la predicción media: \\[\\sum_j \\phi_j(x) = f(x) - E(f(x))\\]\nSi una variable satisface \\(\\phi_W (x) = \\phi_{W\\cup j} (x)\\) para cualquier subconjunto de variables \\(W\\), entonces \\(\\phi_j\\) es cero.\nSi nuestro predictor se puede escribir como \\(f(x) = f_1(x) + f_2(x)\\) entonces la contribución de cada variable en la predicción \\(f(x)\\) es la suma de las contribución en \\(f_1\\) y en \\(f_2\\)\nPara cualquier subconjunto de variables, podemos considerar la contribución de una variable \\(j\\) cuando agregamos al subcojunto la variable \\(j\\) como \\(\\phi_{W \\cup j} (x) - \\phi_W (x)\\). Si para cualquier subconjunto \\(W\\) tenemos que las variables \\(j\\) y \\(k\\) tienen la misma contribución al agregarlos a \\(W\\) (para toda \\(W\\)), entonces \\(\\phi_j =\\phi_k\\).\n\nAhora la pregunta es: ¿cómo definimos las contribuciones para un predictor \\(f\\) más complejo? (por ejemplo, un bosque aleatorio). Resulta que el concepto de contribución o atribución es uno que se estudia en teoría de juegos, y se puede demostrar que hay una sola forma de hacerlo cumpliendo las propiedades señaladas arriba: mediante los valores de Shapley.\nAntes de seguir, veremos las dificultades que podemos encontrar para definir la atribución o contribuciones:\n\nEjemplo\nConsideramos un modelo lineal para empezar. Por ejemplo \\[f(x_1, x_2) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2,\\] ¿Cómo podemos atribuir contribuciones de \\(x_1\\) y \\(x_2\\) para una entrada particular \\(x^* = (x_1^*, x_2^*).\\)\nEn primer lugar, podríamos ver qué pasa cuando pensamos que no tenemos ninguna variable disponible. Esta cantidad la podríamos definir como (promediando sobre la muestra de entrenamiento): \\[v(0) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{(i)}) = \\beta_0 + \\beta_1 \\bar{x_1} + \\beta_2\\bar{x_2}\\]\nque es la predicción media (usamos linealidad de \\(f\\)).\nAhora supongamos que usamos la variable \\(x_1\\). En este caso, podríamos calcular una predicción promediando sobre la variable \\(x_2\\), que no tenemos disponible:\n\\[v(1) = \\frac{1}{N}\\sum_i f(x_1^{*},x_2^{(i)}) = \\beta_0 + \\beta_1 x_1^* + \\beta_2\\bar{x_2}\\] Entonces la contribución cuando agregamos la variable 1 es: \\[v(1) -v(0) = \\beta_1(x_1^* - \\bar{x_1})\\] Ahora supongamos que ya tenemos la variable 2, y agregamos la variable 1. La predicción con todas las variables es \\[v(1,2) = \\beta_0 + \\beta_1 x_1^* + \\beta_2 x_2^*\\] y cuando solo tenemos la variable \\(2\\) es \\[v(2) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{*}) = \\beta_0 + \\beta_1\\bar{x_1} + \\beta_2 x_2^*, \\] y la contribución en este caso es: \\[v(1,2) - v(2) = \\beta_1 (x_1^* - \\bar{x_1}),\\] que es la misma cantidad que \\(v(1)-v(0)\\). Así que no importa como “incluyamos” la variable 1 bajo este criterio, el resultado es el mismo, y podemos definir la contribución de \\(x_1\\) como definimos arriba para un predictor lineal.\nAhora consideremos el caso más complicado donde tenemos una interacción multiplicativa:\n\\[f(x_1, x_2) = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\beta_3x_1x_2\\]\nVeremos qué pasa cuando seguimos la misma idea de arriba: en primer lugar, tenemos\n\\[v(0) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{(i)}) =\n\\beta_0 + \\beta_1 \\overline{x_1} + \\beta_2\\overline{x_2} + \\beta_3\\overline{x_1x_2}\\]\ny \\[v(1) = \\frac{1}{N}\\sum_i f(x_1^{*},x_2^{(i)}) = \\beta_0 + \\beta_1 x_1^* + \\beta_2\\bar{x_2} + (\\beta_3 \\bar{x_2})x_1^*\\] La contribución sería \\[v(1)-v(0) = \\beta_1(x_1^* - \\bar{x_1}) +\\beta_3\\bar{x_2}\\left(x_1^*-\\frac{\\overline{x_1x_2}}{\\bar{x_2}}\\right)\\] Por otro lado, igual que arriba\n\\[v(2) = \\frac{1}{N}\\sum_i f(x_1^{(i)},x_2^{*}) = \\beta_0 + \\beta_1\\bar{x_1} + \\beta_2 x_2^* + \\beta_3x_2^*\\bar{x_1}\\]\ny entonces: \\[v(1,2) - v(2) = \\beta_1(x_1^*-\\bar{x_1})+ \\beta_3x_2^*(x_1^* - \\bar{x_1})\\] Y las cantidades \\(v(1) - v(0)\\) y \\(v(1,2)- v(2)\\) no son iguales en este caso. Tenemos dos órdenes distintos para poner las variables en el modelo, así que ponemos:\n\\[\\phi_1 = \\frac{1}{2}(v(1,2)-v(2)) + \\frac{1}{2}(v(1)-v(0)).\\] y análogamente para \\[\\phi_2 = \\frac{1}{2}(v(1,2)-v(1)) + \\frac{1}{2}(v(2)-v(0)).\\]\nNótese que estas cantidades satisfacen: \\(\\phi_1 + \\phi_2 = v(1,2) -v(0) = f(x_1^*, x_2^*) - \\frac{1}{N}\\sum_i f(x_1^{(i)}, x_2^{(i)}),\\)\nes decir, la diferencia entre la predicción que nos interesa y la predicción media, que es la primera propiedad que buscábamos.\nConsideremos ahora el caso de 3 variables, y nos interesa encontrar la contribución de la variable 1. Ahora tenemos 6 órdenes distintos para poner las variables en el predictor. Dos de ellos terminan con \\(v(1,2,3) - v(2,3)\\), uno de ellos contiene la comparación \\(v(1,2)-v(2)\\), uno de ellos contiene la comparación \\(v(1,3)-v(3)\\) y dos de ellos contiene la comparación \\(v(1)-v(0)\\). Así que ponderaríamos:\n\\[\\phi_1 = \\frac{2}{6}(v(1,2,3) - v(2,3)) + \\frac{1}{6}(v(1,2)-v(2)) + \\frac{1}{6}(v(1,3)-v(3)) + \\frac{2}{6}(v(1)-v(0))\\] Puedes checar que \\(\\phi_1+ \\phi_2 +\\phi_3 = f(x_1^*, x_2^*, x_3^*) - \\frac{1}{N}\\sum_i f(x_1^{(i)}, x_2^{(i)}, x_3^{(i)})\\). Los valores así definidos satisfacen todas las propiedades que discutimos arriba.\n\n\n14.5.1 Valores de Shapley\nLa idea principal en los valores de Shapley es que para entender la contribución de una variable, es necesario considerar su contribución según distintos órdenes en el que consideramos las variables.\nDefinimos primero la contribución de un conjunto de variables \\(S\\subset M,\\) donde \\(M =\\{x_1, x_2\\ldots, x_p\\}\\), marginalizando como sigue con la muestra de entrenamiento:\nSuponiendo que \\(S = \\{1,2,3,\\ldots ,k\\}\\) entonces \\[v(S) = \\frac{1}{N}\\sum_i f(x_1^*,x_2^*,\\ldots, x_k^*, x_k^{(i)}, x_{k+1}^{(i)}, \\ldots, x_p^{(i)})\\] Donde se mantienen fijos los elementos de \\(S\\), y promediamos sobre los elementos que no están en \\(S\\). Esto podemos hacerlo para cualquier subconjunto \\(S\\) de variables. Ahora consideramos una variable dada \\(j\\), y consideramos las diferencias: \\[v(S\\cup {j}) - v(S)\\] Es decir, cuánto contribuye la variable \\(j\\) cuando la agregamos a las variables que ya teníamos en \\(S\\).\nEl valor de Shapley de la variable \\(j\\) es\n\\[\\phi_j = \\sum_{S\\subset M, j\\notin S} \\frac{|S|!(M - |S| -1)!}{M!} (v(S\\cup j) -v(S))\\]\nLos ponderadores están dados como discutimos arriba: hay \\(M!\\) posibles ordenamientos de las variables. Los ordenamientos que incluyen la comparación \\(v(S\\cup j) -v(S)\\) satisfacen: las primeras \\(|S|\\) variables pueden haber entrado de \\(|S|!\\) maneras, luego sigue la \\(j\\), y el restante de las variables pueden entrar de \\((M-|S|-1)!\\) maneras.\n\n\nEjemplo: precios de casas\nConsideremos entonces el bosque que construimos arriba para los precios de casas, y consideramos una casa particular:\n\nset.seed(232)\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[21,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\n\nLa suma de los valores phi da (aproximadamente, por el algoritmo usado) la diferencia entre nuestra predicción para esta instancia y la predicción promedio:\n\nsum(shapley_val$results$phi)\n\n[1] -10.40783\n\n\nAhora hacemos otro ejemplo:\n\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[52,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\nsum(shapley_val$results$phi)\n\n[1] -47.69727\n\n\n\nshapley_val <- Shapley$new(predictor, x.interest = casas_entrena[121,])\nshapley_val$results <- filter(shapley_val$results, feature %in% vars_usadas)\nshapley_val$plot()\n\n\n\n\nObservaciones:\n\nNótese que las contribuciones de cada variable en un mismo nivel puede ser diferente en diferentes instancias, pues los modelos que usamos típicamente incluyen interacciones.\nCalcular sobre todas las posibles permutaciones de las variables es demasiado costoso. Para estimar el valor de Shapley, en iml se toma una muestra de permutaciones, y para cada una se calcula el valor correspondiente \\(v(S\\cup j) - v(S),\\) dependiendo dónde aparece \\(j\\), y promediando sobre las variables que no aparecen en \\(S\\) como mostramos arriba. No es necesario ponderar la muestra de permutaciones.\nExisten también versiones adaptadas a los árboles que son más rápìdas.\n\n\n\nEjemplo\nPara un modelo de clasificación, xgboost calcula los valores de shapley en la escala logit (del predictor lineal)\n\nlibrary(xgboost)\nnse_tbl <- nse_entrena %>% ungroup %>% \n dplyr::select(ab, conex_inte, ocupados, \n num_tvd, educa_jefe, tam_loc, combustible, num_auto)\nmat_x <- model.matrix(ab ~ -1 + ., nse_tbl)\nd_entrena <- xgb.DMatrix(data = mat_x, label = nse_tbl$ab)\nboost_nse <- xgboost(mat_x, label = as.numeric(nse_tbl$ab == \"1\"), \n nrounds = 500, eta = 0.1, \n max_depth = 3, subsample = .5,\n objective = \"binary:logistic\", nthread = 2, verbose = 0)\nshap_nse <- xgb.plot.shap(mat_x, model = boost_nse, top_n = 8, n_col = 2)\n\n\n\npred_base <- mean(nse_tbl$ab==1)\n\n\ndata <- shap_nse$data\ncontrib <- shap_nse$shap_contrib\ndat_g <- tibble(educa_jefe = data[, \"educa_jefe\"], shap_educa = contrib[, \"educa_jefe\"],\n shap_ocupados = contrib[, \"ocupados\"], ) %>% \n sample_n(2000)\nggplot(dat_g, aes(x = educa_jefe, y = shap_educa, colour = shap_ocupados)) + \n geom_jitter(alpha = 0.9, width = 0.5, height = 0) +\n scale_color_gradient(low=\"orange\", high=\"black\")\n\n\n\n\nY podemos ver predicciones individuales:\n\ngraf_caso <- function(data, contrib, ind_caso, pred_base){\n dat_tbl <- tibble(variable = colnames(data), valor = data[ind_caso, ], \n shap = contrib[ind_caso,]) %>% \n unite(\"variable_valor\", variable, valor) %>% \n mutate(variable_valor = fct_reorder(variable_valor, shap))\n pred_logit <- log(pred_base / (1-pred_base)) + sum(dat_tbl$shap)\n pred <- 1 / (1 + exp(-pred_logit))\n sprintf(\"Predicción base: %0.2f, Predicción: %0.2f\", pred_base, pred) %>% print\n ggplot(dat_tbl, aes(x= variable_valor, y = shap, ymax = shap, ymin = 0)) + \n coord_flip() +geom_point() + geom_hline(yintercept = 0, colour = \"gray\") +\n geom_linerange() \n}\ngraf_caso(data, contrib, 10, pred_base)\n\n[1] \"Predicción base: 0.17, Predicción: 0.87\"\n\n\n\n\n\n\ngraf_caso(data, contrib, 102, pred_base)\n\n[1] \"Predicción base: 0.17, Predicción: 0.85\"\n\n\n\n\n\nDiscusión:\n\nUn valor de contribución que puede ser más apropiado para los valores Shapley es condicionando en cada caso a la información que se tiene durante el proceso de adición de variables. Es decir, usamos \\[v(S) = E_{X_C|X_S}\\left [ f(X_S,X_C) \\, | \\, X_S = x_s^* \\right].\\] Esta cantidad es teóricamente más apropiada para hacer predicciones cuando “no tenemos” las variables de \\(X_C\\). Sin embargo, calcular esta cantidad es considerablemente difícil, pues requiere modelar también la conjunta de \\((X_1,\\ldots, X_p)\\), lo cuál en general es difícil. Hasta en regresión lineal, incluso sin interacciones, no es trivial hacer estos cálculos. Generalmente se utilizan como mostramos arriba. Una desventaja clara del proceso como lo mostramos arriba es que puede ser que al hacer estos promedios, usemos partes del modelo con poca información y malas predicciones. Los valores de Shapley pueden ser ruidosos en este caso.\nLos que mostramos arriba son llamados comunmente valores SHAP o explicación aditiva de Shapley. Existen muchas variaciones de la aplicación de valores de Shapley para entender predictores y es un tema de investigación activo.\n\n\n\n\n\nHastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2017. The Elements of Statistical Learning. Springer Series en Statistics. Springer New York Inc. http://web.stanford.edu/~hastie/ElemStatLearn/." }, { "objectID": "81-apendice-descenso.html#cálculo-del-gradiente", @@ -683,14 +683,14 @@ "href": "81-apendice-descenso.html#implementación", "title": "Apéndice A — Apéndice 1: descenso en gradiente", "section": "A.2 Implementación", - "text": "A.2 Implementación\nEn este punto, podemos intentar una implementación simple basada en el código anterior para hacer descenso en gradiente para nuestro problema de regresión (es un buen ejercicio). En lugar de eso, mostraremos cómo usar librerías ahora estándar para hacer esto. En particular usamos keras (con tensorflow), que tienen la ventaja:\n\nEn tensorflow y keras no es necesario calcular las derivadas a mano. Utiliza diferenciación automática, que no es diferenciación numérica ni simbólica: se basa en la regla de la cadena y la codificación explícita de las derivadas de funciones elementales.\n\n\nlibrary(tidymodels)\nlibrary(keras)\nsource(\"../R/casas_traducir_geo.R\")\nset.seed(68821)\n# dividir muestra\ncasas_split <- initial_split(casas |>\n select(precio_m2_miles, area_hab_m2, calidad_gral, num_coches), \n prop = 0.75)\n# obtener muestra de entrenamiento\ncasas_entrena <- training(casas_split)\ncasas_receta <- recipe(precio_m2_miles ~ ., casas_entrena) \n\n\n# definición de estructura del modelo (regresión lineal)\nx_ent <- casas_receta |> prep() |> juice() |> select(-precio_m2_miles) |> as.matrix()\ny_ent <- casas_receta |> prep() |> juice() |> pull(precio_m2_miles)\nn_entrena <- nrow(x_ent)\ncrear_modelo <- function(lr = 0.01){\n modelo_casas <- \n keras_model_sequential() |>\n layer_dense(units = 1, #una sola respuesta,\n activation = \"linear\", # combinar variables linealmente\n kernel_initializer = initializer_constant(0), #inicializamos coeficientes en 0\n bias_initializer = initializer_constant(0)) #inicializamos ordenada en 0\n # compilar seleccionando cantidad a minimizar, optimizador y métricas\n modelo_casas |> compile(\n loss = \"mean_squared_error\", # pérdida cuadrática\n optimizer = optimizer_sgd(learning_rate = lr), # descenso en gradiente\n metrics = list(\"mean_squared_error\"))\n modelo_casas\n}\n# tasa de aprendizaje es lr, tenemos que poner una tasa chica (prueba)\nmodelo_casas <- crear_modelo(lr = 0.00001)\n# Ahora iteramos\n# Primero probamos con un número bajo de iteraciones\nhistoria <- modelo_casas |> fit(\n x_ent, # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 20, # número de iteraciones\n verbose = 0\n)\n\n\nplot(historia, metrics = \"mean_squared_error\", smooth = FALSE) +\n geom_line()\n\n\n\nhistoria$metrics$mean_squared_error |> round(4)\n\n [1] 1.7903 0.7636 0.4549 0.3621 0.3342 0.3258 0.3232 0.3224 0.3222 0.3221\n[11] 0.3220 0.3220 0.3220 0.3219 0.3219 0.3219 0.3218 0.3218 0.3218 0.3217\n\n\nProbamos con más corridas para checar convergencia:\n\n# Agregamos iteraciones: esta historia comienza en los últimos valores de\n# la corrida anterior\nhistoria <- modelo_casas |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 1000, # número de iteraciones\n verbose = 0\n)\n\n\nplot(historia, metrics = \"mean_squared_error\", smooth = FALSE) \n\n\n\n\nEl modelo parece todavía ir mejorando. Veamos de todas formas los coeficientes estimados hasta ahora:\n\nkeras::get_weights(modelo_casas)\n\n[[1]]\n [,1]\n[1,] 0.007331367\n[2,] 0.016437242\n[3,] 0.004633877\n\n[[2]]\n[1] 0.00302865\n\n\nLa implementación oficial de R es lm, que en general tiene buen desempeño para datos que caben en memoria:\n\nlm(precio_m2_miles ~ area_hab_m2 + calidad_gral + num_coches, \n data = casas_entrena) |> \n coef()\n\n (Intercept) area_hab_m2 calidad_gral num_coches \n 0.66869194 -0.00449751 0.16807663 0.13115749 \n\n\nDe modo que todavía requerimos más iteraciones para alcanzar convergencia. ¿Por qué la convergencia es tan lenta? En parte, la razón es que las escalas de las variables de entrada son muy diferentes, de modo que es difícil ajustar una tasa de aprendizaje constante que funcione bien. Podemos remediar esto poniendo todas las entradas en la misma escala (normalizando)" + "text": "A.2 Implementación\nEn este punto, podemos intentar una implementación simple basada en el código anterior para hacer descenso en gradiente para nuestro problema de regresión (es un buen ejercicio). En lugar de eso, mostraremos cómo usar librerías ahora estándar para hacer esto. En particular usamos keras (con tensorflow), que tienen la ventaja:\n\nEn tensorflow y keras no es necesario calcular las derivadas a mano. Utiliza diferenciación automática, que no es diferenciación numérica ni simbólica: se basa en la regla de la cadena y la codificación explícita de las derivadas de funciones elementales.\n\n\nlibrary(tidymodels)\nlibrary(keras)\nsource(\"../R/casas_traducir_geo.R\")\nset.seed(68821)\n# dividir muestra\ncasas_split <- initial_split(casas |>\n select(precio_m2_miles, area_hab_m2, calidad_gral, num_coches), \n prop = 0.75)\n# obtener muestra de entrenamiento\ncasas_entrena <- training(casas_split)\ncasas_receta <- recipe(precio_m2_miles ~ ., casas_entrena) \n\n\n# definición de estructura del modelo (regresión lineal)\nx_ent <- casas_receta |> prep() |> juice() |> select(-precio_m2_miles) |> as.matrix()\ny_ent <- casas_receta |> prep() |> juice() |> pull(precio_m2_miles)\nn_entrena <- nrow(x_ent)\ncrear_modelo <- function(lr = 0.01){\n modelo_casas <- \n keras_model_sequential() |>\n layer_dense(units = 1, #una sola respuesta,\n activation = \"linear\", # combinar variables linealmente\n kernel_initializer = initializer_constant(0), #inicializamos coeficientes en 0\n bias_initializer = initializer_constant(0)) #inicializamos ordenada en 0\n # compilar seleccionando cantidad a minimizar, optimizador y métricas\n modelo_casas |> compile(\n loss = \"mean_squared_error\", # pérdida cuadrática\n optimizer = optimizer_sgd(learning_rate = lr), # descenso en gradiente\n metrics = list(\"mean_squared_error\"))\n modelo_casas\n}\n# tasa de aprendizaje es lr, tenemos que poner una tasa chica (prueba)\nmodelo_casas <- crear_modelo(lr = 0.00001)\n# Ahora iteramos\n# Primero probamos con un número bajo de iteraciones\nhistoria <- modelo_casas |> fit(\n x_ent, # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 20, # número de iteraciones\n verbose = 0\n)\n\n\nplot(historia, metrics = \"mean_squared_error\", smooth = FALSE) +\n geom_line()\n\n\n\nhistoria$metrics$mean_squared_error |> round(4)\n\n [1] 1.7903 0.7636 0.4549 0.3621 0.3342 0.3258 0.3232 0.3224 0.3222 0.3221\n[11] 0.3220 0.3220 0.3220 0.3219 0.3219 0.3219 0.3218 0.3218 0.3218 0.3217\n\n\nProbamos con más corridas para checar convergencia:\n\n# Agregamos iteraciones: esta historia comienza en los últimos valores de\n# la corrida anterior\nhistoria <- modelo_casas |> fit(\n as.matrix(x_ent), # x entradas\n y_ent, # y salida o target\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = 1000, # número de iteraciones\n verbose = 0\n)\n\n\nplot(historia, metrics = \"mean_squared_error\", smooth = FALSE) \n\n\n\n\nEl modelo parece todavía ir mejorando. Veamos de todas formas los coeficientes estimados hasta ahora:\n\nkeras::get_weights(modelo_casas)\n\n[[1]]\n [,1]\n[1,] 0.007331367\n[2,] 0.016437238\n[3,] 0.004633876\n\n[[2]]\n[1] 0.003028649\n\n\nLa implementación oficial de R es lm, que en general tiene buen desempeño para datos que caben en memoria:\n\nlm(precio_m2_miles ~ area_hab_m2 + calidad_gral + num_coches, \n data = casas_entrena) |> \n coef()\n\n (Intercept) area_hab_m2 calidad_gral num_coches \n 0.66869194 -0.00449751 0.16807663 0.13115749 \n\n\nDe modo que todavía requerimos más iteraciones para alcanzar convergencia. ¿Por qué la convergencia es tan lenta? En parte, la razón es que las escalas de las variables de entrada son muy diferentes, de modo que es difícil ajustar una tasa de aprendizaje constante que funcione bien. Podemos remediar esto poniendo todas las entradas en la misma escala (normalizando)" }, { "objectID": "81-apendice-descenso.html#normalización-de-entradas", "href": "81-apendice-descenso.html#normalización-de-entradas", "title": "Apéndice A — Apéndice 1: descenso en gradiente", "section": "A.3 Normalización de entradas", - "text": "A.3 Normalización de entradas\nLa convergencia de descenso en gradiente (y también el desempeño numérico para otros algoritmos) puede dificultarse cuando las variables tienen escalas muy diferentes. Esto produce curvaturas altas en la función que queremos minimizar.\nEn este ejemplo simple, una variable tiene desviación estándar 10 y otra 1:\n\nx1 <- rnorm(100, 0, 5) \nx2 <- rnorm(100, 0, 1) + 0.1*x1\ny <- 0*x1 + 0*x2 + rnorm(100, 0, 0.1) \ndat <- tibble(x1, x2, y)\nrss <- function(beta) mean((as.matrix(dat[, 1:2]) %*% beta - y)^2) \ngrid_beta <- expand.grid(beta1 = seq(-1, 1, length.out = 50), \n beta2 = seq(-1, 1, length.out = 50))\nrss_1 <- apply(grid_beta, 1, rss) \ndat_x <- data.frame(grid_beta, rss_1)\nggplot(dat_x, aes(x = beta1, y = beta2, z = rss_1)) + \n geom_contour(binwidth = 0.5) +\n coord_equal() \n\n\n\n\nEn algunas direcciones el gradiente es muy grande, y en otras chico. Esto implica que la convergencia puede ser muy lenta en algunas direcciones, puede diverger en otras, y que hay que ajustar el paso \\(\\eta > 0\\) con cuidado, dependiendo de dónde comiencen las iteraciones.\nPor ejemplo, con un tamaño de paso relativamente chico, damos unos saltos grandes al principio y luego avanzamos muy lentamente:\n\ngrad_calc <- function(x_ent, y_ent){\n # calculamos directamente el gradiente\n salida_grad <- function(beta){\n n <- length(y_ent)\n f_beta <- as.matrix(cbind(1, x_ent)) %*% beta\n e <- y_ent - f_beta\n grad_out <- - as.numeric(t(cbind(1, x_ent)) %*% e) / n\n names(grad_out) <- c('Intercept', colnames(x_ent))\n grad_out\n }\n salida_grad\n}\ngrad_sin_norm <- grad_calc(dat[, 1:2, drop = FALSE], dat$y)\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.02, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\nSi incrementamos el tamaño de paso observamos también convergencia lenta. En este caso particular, subir más el tamaño de paso puede producir divergencia:\n\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.07, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\nUna normalización usual es con la media y desviación estándar, donde hacemos, para cada variable de entrada \\(j=1,2,\\ldots, p\\) \\[ x_j^{(i)} = \\frac{ x_j^{(i)} - \\bar{x}_j}{s_j}\\] donde \\[\\bar{x}_j = \\frac{1}{N} \\sum_{i=1}^N x_j^{(i)}\\] \\[s_j = \\sqrt{\\frac{1}{N-1}\\sum_{i=1}^N (x_j^{(i)}- \\bar{x}_j )^2}\\] es decir, centramos y normalizamos por columna. Otra opción común es restar el mínimo y dividir entre la diferencia del máximo y el mínimo, de modo que las variables resultantes toman valores en \\([0,1]\\).\nEntonces escalamos antes de ajustar:\n\nx1_s = (x1 - mean(x1))/sd(x1)\nx2_s = (x2 - mean(x2))/sd(x2)\ndat <- tibble(x1_s = x1_s, x2_s = x2_s, y = y)\nrss <- function(beta) mean((as.matrix(dat[, 1:2]) %*% beta - y)^2) \ngrid_beta <- expand.grid(beta1 = seq(-1, 1, length.out = 50), \n beta2 = seq(-1, 1, length.out = 50))\nrss_1 <- apply(grid_beta, 1, rss) \ndat_x <- data.frame(grid_beta, rss_1)\nggplot(dat_x, aes(x = beta1, y = beta2, z = rss_1)) + \n geom_contour(binwidth = 0.5) +\n coord_equal() \n\n\n\n\nNótese que los coeficientes ajustados serán diferentes a los del caso no normalizado.\nSi normalizamos, obtenemos convergencia más rápida\n\ngrad_sin_norm <- grad_calc(dat[, 1:2, drop = FALSE], dat$y)\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.5, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\n\n\n\n\n\n\nTip\n\n\n\nCuando normalizamos antes de ajustar el modelo, las predicciones deben hacerse con entradas normalizadas. La normalización se hace con los mismos valores que se usaron en el entrenamiento (y no recalculando medias y desviaciones estándar con el conjunto de prueba). En cuanto a la forma funcional del predictor \\(f\\), el problema con entradas normalizadas es equivalente al de las entradas no normalizadas. Asegúrate de esto escribiendo cómo correponden los coeficientes de cada modelo normalizado con los coeficientes del modelo no normalizado.\n\n\nSupongamos que el modelo en las variables originales es \\[{f}_\\beta (X) = \\beta_0 + \\beta_1 X_1 + \\beta_2 X_2 + \\cdots + \\beta_p X_p,\\] Consideramos el modelo con variables estandarizadas \\[{g}_{\\beta^s} (X) = \\beta_0^s + \\beta_1^s Z_1 + \\beta_2^s Z_2 + \\cdots + \\beta_p^s Z_p,\\]\nSustituyendo \\(Z_j = (X_j - \\mu_j)/s_j,\\)\n\\[{g}_{\\beta^s} (X) = (\\beta_0^s - \\sum_{j=1}^p \\beta_j^s \\mu_j/s_j) + \\frac{\\beta_1^s}{s_j} X_1 + \\frac{\\beta_2^s}{s_2} X_2 + \\cdots + \\frac{\\beta_p^s}{s_p} X_p,\\] Y vemos que tiene la misma forma funcional de \\(f_\\beta(X)\\). Si la solución de mínimos cuadrados es única, entonces una vez que ajustemos tenemos que tener \\(\\hat{f}_\\beta(X) = \\hat{g}_{\\beta^s} (X)\\), lo que implica que \\[\\hat{\\beta}_0 = \\hat{\\beta}_0^s - \\sum_{j=1}^p \\hat{\\beta}_j^s\\mu_j/s_j\\] y \\[\\hat{\\beta}_j = \\hat{\\beta}_j^s/s_j.\\]\nNótese que para pasar del problema estandarizado al no estandarizado simplemente se requiere escalar los coeficientes por la \\(s_j\\) correspondiente.\n\nEjemplo\nRepetimos nuestro modelo, pero normalizando las entradas:\n\n# usamos recipes para este ejemplo, no necesitas usarlo\ncasas_receta <- recipe(precio_m2_miles ~ ., casas_entrena) |>\n step_normalize(all_predictors()) \ncasas_receta |> summary()\n\n# A tibble: 4 × 4\n variable type role source \n <chr> <list> <chr> <chr> \n1 area_hab_m2 <chr [2]> predictor original\n2 calidad_gral <chr [2]> predictor original\n3 num_coches <chr [2]> predictor original\n4 precio_m2_miles <chr [2]> outcome original\n\n\n\nmodelo_lineal <- linear_reg() |>\n set_engine(\"lm\")\ncasas_flujo <- workflow() |>\n add_recipe(casas_receta) |> \n add_model(modelo_lineal)\n\n\nlibrary(keras)\n# definición de estructura del modelo (regresión lineal)\nx_ent_s <- prep(casas_receta) |> juice() |> select(-precio_m2_miles) |> \n as.matrix()\najustar_casas <- function(modelo, x, y, n_epochs = 100){\n ajuste <- modelo |> fit(\n as.matrix(x), y,\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = n_epochs, # número de iteraciones\n verbose = 0) |> as_tibble()\n ajuste\n}\nmodelo_casas_ns <- crear_modelo(0.00001)\nmodelo_casas_s <- crear_modelo(0.2)\nhistoria_s <- ajustar_casas(modelo_casas_s, x_ent_s, y_ent) |>\n mutate(tipo = \"Estandarizar\")\nhistoria_ns <- ajustar_casas(modelo_casas_ns, x_ent, y_ent) |> \n mutate(tipo = \"Sin estandarizar\")\nhistoria <- bind_rows(historia_ns, historia_s) |> filter(metric == \"mean_squared_error\")\nggplot(historia, aes(x = epoch, y = value, colour = tipo)) +\n geom_line() + geom_point() +scale_x_log10() + scale_y_log10()\n\n\n\n\nObservamos que el modelo con datos estandarizados convergió:\n\nkeras::get_weights(modelo_casas_s)\n\n[[1]]\n [,1]\n[1,] -0.22045916\n[2,] 0.23149671\n[3,] 0.09673478\n\n[[2]]\n[1] 1.295027\n\ncoef(lm.fit(cbind(1,x_ent_s), y_ent))\n\n area_hab_m2 calidad_gral num_coches \n 1.29502679 -0.22045919 0.23149675 0.09673476 \n\n\nMientras que el modelo no estandarizado todavía requiere iteraciones:\n\nkeras::get_weights(modelo_casas_ns)\n\n[[1]]\n [,1]\n[1,] 0.0079817399\n[2,] 0.0019486139\n[3,] 0.0005532746\n\n[[2]]\n[1] 0.0003491883\n\ncoef(lm.fit(cbind(1, x_ent), y_ent))\n\n area_hab_m2 calidad_gral num_coches \n 0.66869194 -0.00449751 0.16807663 0.13115749" + "text": "A.3 Normalización de entradas\nLa convergencia de descenso en gradiente (y también el desempeño numérico para otros algoritmos) puede dificultarse cuando las variables tienen escalas muy diferentes. Esto produce curvaturas altas en la función que queremos minimizar.\nEn este ejemplo simple, una variable tiene desviación estándar 10 y otra 1:\n\nx1 <- rnorm(100, 0, 5) \nx2 <- rnorm(100, 0, 1) + 0.1*x1\ny <- 0*x1 + 0*x2 + rnorm(100, 0, 0.1) \ndat <- tibble(x1, x2, y)\nrss <- function(beta) mean((as.matrix(dat[, 1:2]) %*% beta - y)^2) \ngrid_beta <- expand.grid(beta1 = seq(-1, 1, length.out = 50), \n beta2 = seq(-1, 1, length.out = 50))\nrss_1 <- apply(grid_beta, 1, rss) \ndat_x <- data.frame(grid_beta, rss_1)\nggplot(dat_x, aes(x = beta1, y = beta2, z = rss_1)) + \n geom_contour(binwidth = 0.5) +\n coord_equal() \n\n\n\n\nEn algunas direcciones el gradiente es muy grande, y en otras chico. Esto implica que la convergencia puede ser muy lenta en algunas direcciones, puede diverger en otras, y que hay que ajustar el paso \\(\\eta > 0\\) con cuidado, dependiendo de dónde comiencen las iteraciones.\nPor ejemplo, con un tamaño de paso relativamente chico, damos unos saltos grandes al principio y luego avanzamos muy lentamente:\n\ngrad_calc <- function(x_ent, y_ent){\n # calculamos directamente el gradiente\n salida_grad <- function(beta){\n n <- length(y_ent)\n f_beta <- as.matrix(cbind(1, x_ent)) %*% beta\n e <- y_ent - f_beta\n grad_out <- - as.numeric(t(cbind(1, x_ent)) %*% e) / n\n names(grad_out) <- c('Intercept', colnames(x_ent))\n grad_out\n }\n salida_grad\n}\ngrad_sin_norm <- grad_calc(dat[, 1:2, drop = FALSE], dat$y)\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.02, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\nSi incrementamos el tamaño de paso observamos también convergencia lenta. En este caso particular, subir más el tamaño de paso puede producir divergencia:\n\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.07, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\nUna normalización usual es con la media y desviación estándar, donde hacemos, para cada variable de entrada \\(j=1,2,\\ldots, p\\) \\[ x_j^{(i)} = \\frac{ x_j^{(i)} - \\bar{x}_j}{s_j}\\] donde \\[\\bar{x}_j = \\frac{1}{N} \\sum_{i=1}^N x_j^{(i)}\\] \\[s_j = \\sqrt{\\frac{1}{N-1}\\sum_{i=1}^N (x_j^{(i)}- \\bar{x}_j )^2}\\] es decir, centramos y normalizamos por columna. Otra opción común es restar el mínimo y dividir entre la diferencia del máximo y el mínimo, de modo que las variables resultantes toman valores en \\([0,1]\\).\nEntonces escalamos antes de ajustar:\n\nx1_s = (x1 - mean(x1))/sd(x1)\nx2_s = (x2 - mean(x2))/sd(x2)\ndat <- tibble(x1_s = x1_s, x2_s = x2_s, y = y)\nrss <- function(beta) mean((as.matrix(dat[, 1:2]) %*% beta - y)^2) \ngrid_beta <- expand.grid(beta1 = seq(-1, 1, length.out = 50), \n beta2 = seq(-1, 1, length.out = 50))\nrss_1 <- apply(grid_beta, 1, rss) \ndat_x <- data.frame(grid_beta, rss_1)\nggplot(dat_x, aes(x = beta1, y = beta2, z = rss_1)) + \n geom_contour(binwidth = 0.5) +\n coord_equal() \n\n\n\n\nNótese que los coeficientes ajustados serán diferentes a los del caso no normalizado.\nSi normalizamos, obtenemos convergencia más rápida\n\ngrad_sin_norm <- grad_calc(dat[, 1:2, drop = FALSE], dat$y)\niteraciones <- descenso(10, c(0, -0.25, -0.75), 0.5, grad_sin_norm)\nggplot(dat_x) + \n geom_contour(aes(x = beta1, y = beta2, z = rss_1), binwidth = 0.5) +\n coord_equal() +\n geom_path(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red') +\n geom_point(data = data.frame(iteraciones[, 2:3]), aes(x=X1, y=X2), colour = 'red')\n\n\n\n\n\n\n\n\n\n\nTip\n\n\n\nCuando normalizamos antes de ajustar el modelo, las predicciones deben hacerse con entradas normalizadas. La normalización se hace con los mismos valores que se usaron en el entrenamiento (y no recalculando medias y desviaciones estándar con el conjunto de prueba). En cuanto a la forma funcional del predictor \\(f\\), el problema con entradas normalizadas es equivalente al de las entradas no normalizadas. Asegúrate de esto escribiendo cómo correponden los coeficientes de cada modelo normalizado con los coeficientes del modelo no normalizado.\n\n\nSupongamos que el modelo en las variables originales es \\[{f}_\\beta (X) = \\beta_0 + \\beta_1 X_1 + \\beta_2 X_2 + \\cdots + \\beta_p X_p,\\] Consideramos el modelo con variables estandarizadas \\[{g}_{\\beta^s} (X) = \\beta_0^s + \\beta_1^s Z_1 + \\beta_2^s Z_2 + \\cdots + \\beta_p^s Z_p,\\]\nSustituyendo \\(Z_j = (X_j - \\mu_j)/s_j,\\)\n\\[{g}_{\\beta^s} (X) = (\\beta_0^s - \\sum_{j=1}^p \\beta_j^s \\mu_j/s_j) + \\frac{\\beta_1^s}{s_j} X_1 + \\frac{\\beta_2^s}{s_2} X_2 + \\cdots + \\frac{\\beta_p^s}{s_p} X_p,\\] Y vemos que tiene la misma forma funcional de \\(f_\\beta(X)\\). Si la solución de mínimos cuadrados es única, entonces una vez que ajustemos tenemos que tener \\(\\hat{f}_\\beta(X) = \\hat{g}_{\\beta^s} (X)\\), lo que implica que \\[\\hat{\\beta}_0 = \\hat{\\beta}_0^s - \\sum_{j=1}^p \\hat{\\beta}_j^s\\mu_j/s_j\\] y \\[\\hat{\\beta}_j = \\hat{\\beta}_j^s/s_j.\\]\nNótese que para pasar del problema estandarizado al no estandarizado simplemente se requiere escalar los coeficientes por la \\(s_j\\) correspondiente.\n\nEjemplo\nRepetimos nuestro modelo, pero normalizando las entradas:\n\n# usamos recipes para este ejemplo, no necesitas usarlo\ncasas_receta <- recipe(precio_m2_miles ~ ., casas_entrena) |>\n step_normalize(all_predictors()) \ncasas_receta |> summary()\n\n# A tibble: 4 × 4\n variable type role source \n <chr> <list> <chr> <chr> \n1 area_hab_m2 <chr [2]> predictor original\n2 calidad_gral <chr [2]> predictor original\n3 num_coches <chr [2]> predictor original\n4 precio_m2_miles <chr [2]> outcome original\n\n\n\nmodelo_lineal <- linear_reg() |>\n set_engine(\"lm\")\ncasas_flujo <- workflow() |>\n add_recipe(casas_receta) |> \n add_model(modelo_lineal)\n\n\nlibrary(keras)\n# definición de estructura del modelo (regresión lineal)\nx_ent_s <- prep(casas_receta) |> juice() |> select(-precio_m2_miles) |> \n as.matrix()\najustar_casas <- function(modelo, x, y, n_epochs = 100){\n ajuste <- modelo |> fit(\n as.matrix(x), y,\n batch_size = nrow(x_ent), # para descenso en gradiente\n epochs = n_epochs, # número de iteraciones\n verbose = 0) |> as_tibble()\n ajuste\n}\nmodelo_casas_ns <- crear_modelo(0.00001)\nmodelo_casas_s <- crear_modelo(0.2)\nhistoria_s <- ajustar_casas(modelo_casas_s, x_ent_s, y_ent) |>\n mutate(tipo = \"Estandarizar\")\nhistoria_ns <- ajustar_casas(modelo_casas_ns, x_ent, y_ent) |> \n mutate(tipo = \"Sin estandarizar\")\nhistoria <- bind_rows(historia_ns, historia_s) |> filter(metric == \"mean_squared_error\")\nggplot(historia, aes(x = epoch, y = value, colour = tipo)) +\n geom_line() + geom_point() +scale_x_log10() + scale_y_log10()\n\n\n\n\nObservamos que el modelo con datos estandarizados convergió:\n\nkeras::get_weights(modelo_casas_s)\n\n[[1]]\n [,1]\n[1,] -0.22045918\n[2,] 0.23149671\n[3,] 0.09673478\n\n[[2]]\n[1] 1.295027\n\ncoef(lm.fit(cbind(1,x_ent_s), y_ent))\n\n area_hab_m2 calidad_gral num_coches \n 1.29502679 -0.22045919 0.23149675 0.09673476 \n\n\nMientras que el modelo no estandarizado todavía requiere iteraciones:\n\nkeras::get_weights(modelo_casas_ns)\n\n[[1]]\n [,1]\n[1,] 0.0079817399\n[2,] 0.0019486138\n[3,] 0.0005532746\n\n[[2]]\n[1] 0.0003491882\n\ncoef(lm.fit(cbind(1, x_ent), y_ent))\n\n area_hab_m2 calidad_gral num_coches \n 0.66869194 -0.00449751 0.16807663 0.13115749" }, { "objectID": "82-apendice-descenso-estocastico.html#algoritmo-de-descenso-estocástico", @@ -725,7 +725,7 @@ "href": "82-apendice-descenso-estocastico.html#ajuste-de-redes-con-descenso-estocástico", "title": "Apéndice B — Apéndice 2: Descenso estocástico", "section": "B.5 Ajuste de redes con descenso estocástico", - "text": "B.5 Ajuste de redes con descenso estocástico\n\nlibrary(keras)\n\n\nset.seed(21321)\nx_ent <- as.matrix(dat_ent[,c('x_1','x_2','x_3')])\nx_valid <- as.matrix(dat_valid[,c('x_1','x_2','x_3')])\ny_ent <- dat_ent$y\ny_valid <- dat_valid$y\n\nEmpezamos con regresión (sin capas ocultas), que se escribe y ajusta como sigue:\n\nmodelo <- keras_model_sequential() \nmodelo |>\n layer_dense(units = 1, \n activation = \"linear\",\n input_shape = c(3))\n\nmodelo |> compile(loss = 'mse',\n optimizer = optimizer_sgd(learning_rate = 0.1, momentum = 0,\n decay = 0))\n\nhistory <- modelo |> \n fit(x_ent, y_ent, \n epochs = 50, batch_size = 10, \n verbose = 0,\n validation_data = list(x_valid, y_valid))\n\nPodemos ver el progreso del algoritmo por época\n\naprendizaje <- as_tibble(history)\nggplot(aprendizaje, \n aes(x=epoch, y=value, colour=data, group=data)) +\n facet_wrap(~metric, ncol = 1) + geom_line() + geom_point(size = 0.5)\n\n\n\n\nVer los pesos:\n\nget_weights(modelo)\n\n[[1]]\n [,1]\n[1,] -2.515035629\n[2,] -0.010236964\n[3,] -0.007203732\n\n[[2]]\n[1] -0.3029693\n\n\nY verificamos que concuerda con la salida de lm:\n\nmod_lineal <- lm(y ~ x_1 + x_2+ x_3, data = dat_ent) \ncoef(mod_lineal)\n\n(Intercept) x_1 x_2 x_3 \n-0.36904266 -2.46877687 -0.07368414 0.06632769 \n\n\n\n\n\n\nGoodfellow, Ian, Yoshua Bengio, y Aaron Courville. 2016. Deep Learning. MIT Press." + "text": "B.5 Ajuste de redes con descenso estocástico\n\nlibrary(keras)\n\n\nset.seed(21321)\nx_ent <- as.matrix(dat_ent[,c('x_1','x_2','x_3')])\nx_valid <- as.matrix(dat_valid[,c('x_1','x_2','x_3')])\ny_ent <- dat_ent$y\ny_valid <- dat_valid$y\n\nEmpezamos con regresión (sin capas ocultas), que se escribe y ajusta como sigue:\n\nmodelo <- keras_model_sequential() \nmodelo |>\n layer_dense(units = 1, \n activation = \"linear\",\n input_shape = c(3))\n\nmodelo |> compile(loss = 'mse',\n optimizer = optimizer_sgd(learning_rate = 0.1, momentum = 0,\n decay = 0))\n\nhistory <- modelo |> \n fit(x_ent, y_ent, \n epochs = 50, batch_size = 10, \n verbose = 0,\n validation_data = list(x_valid, y_valid))\n\nPodemos ver el progreso del algoritmo por época\n\naprendizaje <- as_tibble(history)\nggplot(aprendizaje, \n aes(x=epoch, y=value, colour=data, group=data)) +\n facet_wrap(~metric, ncol = 1) + geom_line() + geom_point(size = 0.5)\n\n\n\n\nVer los pesos:\n\nget_weights(modelo)\n\n[[1]]\n [,1]\n[1,] -2.64986229\n[2,] -0.05415167\n[3,] 0.04267785\n\n[[2]]\n[1] -0.4793204\n\n\nY verificamos que concuerda con la salida de lm:\n\nmod_lineal <- lm(y ~ x_1 + x_2+ x_3, data = dat_ent) \ncoef(mod_lineal)\n\n(Intercept) x_1 x_2 x_3 \n-0.36904266 -2.46877687 -0.07368414 0.06632769 \n\n\n\n\n\n\nGoodfellow, Ian, Yoshua Bengio, y Aaron Courville. 2016. Deep Learning. MIT Press." }, { "objectID": "99-referencias.html",