Validación Cruzada

A veces hay en el mundo una variable aleatoria $X$ que distintos agentes tratan de predecir. Pero en otras ocasiones hay una variable aleatoria $X$, ciertos agentes crean predicciones $Y$ de $X$ y otros están interesados en la predicción de $Y$, tal como describe Matt Levine en Weather Prediction Prediction. The Economist describe los esfuerzos de la India por recuperar la credibilidad de sus estadísticas nacionales tras años de interferencia política. ...

I. Juan Cambeiro escribe en Asterisk What Comes After COVID. El covid nos aburre y no nos interesa, pero el artículo es un ejercicio de “probabilidad aplicada” —en el que se estudia cuándo y qué causará la próxima pandemia, pero eso es casi lo de menos— del que muchos podrán sacar provecho. II. La mayor parte de los artículos en economía son inútiles; todos los involucrados lo saben. Fuera del primer cuartil, todo es esencialmente una estafa que no sobreviviría una revisión crítica." ...

He escrito esta entrada como una introducción a lo que se cuenta aquí, aquí y aquí sobre el asunto de la relación entre la optimización (como parte del proceso de ajuste de modelos) y la generalización (o su capacidad para aprender sobre el mundo y no solo sobre los datos de entrenamiento). En los enlaces, el lector encontrará planteadas una serie de cuestiones sobre cómo y por qué generalizan los (o cierto tipo de) modelos en lugar de, simplemente, no hacerlo. ...

Supongamos que tenemos un modelo construido sobre unos datos $(x_i, y_i)$. Para cada $x_i$, el valor $y_i$ es una realización de una variable aleatoria $Y_i$ con distribución $F_i(y)$. Por simplificar, podemos suponer, además, que para el ajuste se utiliza el error cuadrático. Entonces, lo mejor que puede hacer el modelo es encontrar la media $\mu_i$ de cada $Y_i$ —bueno, en realidad, querría encontrar $\mu_x$ para cada $x$ potencial, pero hoy vamos a dejar esa discusión aparcada—. ...

Al construir modelos, queremos minimizar $$ l(\theta) = \int L(y, f_\theta(x)) dP(x,y),$$ donde $L$ es una determinada función de pérdida (y no, no me refiero exclusivamente a la que tiene un numerillo 2). Pero como de $P(x,y)$ solo conocemos una muestra $(x_i, y_i)$ (dejadme aprovechar la ocasión para utilizar una de mis palabras favoritas: $P(x,y)$ es incognoscible), hacemos uso de la aproximación $$ \int f(x) dP(x) \approx \frac{1}{N} \sum f(x_i)$$ ...

Es una pregunta legítima —en el sentido de que ignoro la respuesta— que tengo. Para plantearla en sus debidos términos: Contexto: Tenemos modelos y queremos compararlos. Queremos que funcionen en el universo, pero solo disponemos de él una muestra. Acto 1: Para desatascar el nudo lógico, recurrimos a técnicas como: Entrenamiento y validación, jackknife y sobre todo, su popular evolución, la validación cruzada. Todas ellas bien sabidas y discutidas en todos los manuales. ...

Está de moda usar caret para estas cosas, pero yo estoy todavía acostumbrado a hacerlas a mano. Creo, además, que es poco instructivo ocultar estas cuestiones detrás de funciones de tipo caja-negra-maravillosa a quienes se inician en el mundo de la construcción y comparación de modelos. Muestro, por tanto, código bastante simple para la validación cruzada de un modelo con R: # genero ids ids <- rep(1:10, length.out = nrow(cars)) # Nota: da igual si nrow(df) no es múltiplo de 10 # los aleatorizo ids <- sample(ids) # esto devuelve una lista de dfs: preds.cv <- lapply(unique(ids), function(i){ preds <- predict(lm(dist ~ speed, data = cars[ids != i,]), cars[ids == i,]) data.frame( preds = preds, real = cars[ids == i,]$dist) }) # "apilo" los dfs: preds.cv <- do.call(rbind, preds.cv) # calculo el rmse rmse <- sqrt(mean((preds.cv$preds - preds.cv$real)^2)) Sí, estoy usando el RMSE aunque sea un detractor del mismo. ...

Estoy sin tiempo, así que os suelto el código y me largo a casa a no cenar. Es así: library(parallel) cl <- makeCluster(8) # solo si hay aleatorización # clusterSetRNGStream(cl, 123) clusterEvalQ(cl, { # las librerías necesarias tienen que cargarse # en cada esclavo library(rpart) # en la práctica, hay que cargar los datos # (¿desde fichero?) en cada esclavo my.data <- iris # lo mismo con las funciones necesarias foo <- function(x, dat){ train <- 1:nrow(dat) %% 10 != 1 mod <- rpart(Species ~ ., data = dat[train,]) res <- predict(mod, dat[!train,]) } }) res <- parSapply(cl, 0:9, function(x) foo(x, my.data), simplify = F)