Validación Cruzada

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 es 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 $latex P(x,y)$ solo conocemos una muestra $latex (x_i, y_i)$ (dejadme aprovechar la ocasión para utilizar una de mis palabras favoritas: $latex 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,j
• 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))

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)