Estadística

Me piden bibliografía para unos cursos de ciencia de datos. En particular, de estadística básica. Un texto que reúna los conceptos fundamentales de la cosa para quienes o no los aprendieron en su día o los olvidaron por el camino. Tiene que cumplir algunos requisitos mínimos:

• Que presente los gráficos estadísticos básicos y que no estén construidos con Excel (en 3D).
• Que, a lo más, incluya un único gráfico de tarta.
• Que no sea muy pesado matemáticamente.
• Que sea breve, pero no demasiado.
• Que esté accesible, idealmente en internet, gratuita y legalmente.

Finalmente, si está escrito escrito en español y usa R, mejor aún.

Voy a constuir unos datos artificiales y un modelo de clasificación binaria,

library(mgcv)
library(ggplot2)
library(pROC)

n <- 10000
dat <- gamSim(1, n=n, dist="binary", scale=.33)

lr.fit <- gam(y ~ s(x0, bs="cr") +
    s(x1, bs="cr") + s(x2, bs="cr") +
    s(x3, bs="cr"),
    family=binomial, data=dat,
    method="REML")

y luego (mal hecho: debería hacerlo sobre un conjunto de validación distinto) a obtener las predicciones para las observaciones

res <- data.frame(real = factor(dat$y),
    prob = predict(lr.fit, type = "response"))

que

ggplot(res, aes(x=prob, fill=real)) +
    geom_density(alpha=.3)

representa así:

Me pregunto si el clasificador construido es bueno. Para lo cual voy a construir la curva ROC con

Sí, de drogas de las que mantienen despierto al lumpenazgo. Porque he encontrado (aquí) un conjunto datos muy interesante sobre la valoración que una serie de personas, unas 900, da a una serie de drogas más o menos legales que se llaman —me acabo de enterar— nootrópicos.

El gráfico

extraído de la página enlazada más arriba resume parte de los resultados. No obstante, es sabido entre los que se dedican a los sistemas de recomendación que hay usuarios que tienden a valorar sistemáticamente por encima de la media y otros, por debajo. En los manuales de la cosa suelen recogerse mecanismos más o menos sofisticados para mitigar ese efecto y normalizar las valoraciones entre usuarios. Generalmente, solo exigen matemáticas de bachillerato. Y son meras aproximaciones que no tienen en cuenta circunstancias tales como que puede que un usuario da valoraciones bajas solo porque evalúa productos malos, etc.

Estos días he estado haciendo de campaña promoviendo el uso de nuevas técnicas de análisis de datos en ámbitos como, p.e., el riesgo de crédito, uno de esos campos sujetos al parecer de un regulador (el Banco de España, en este caso).

La gente con la que he debatido al respecto tiende a aplicar esa forma cuasiperfecta de censura que es la autocensura previa. La autocensura previa ni siquiera requiere la acción explícita del censor: es el potencial censurado el que la aplica de mejor o peor gana automáticamente… por si las moscas.

Imaginemos un banco que construye modelos para determinar si se concede o no un crédito. Este banco tiene varias opciones para crear el modelo. Sin embargo, en algunos países el regulador exige que el banco pueda explicar el motivo de la denegación de un crédito cuando un cliente lo solicite.

Esa restricción impediría potencialmente usar modelos de caja negra como el que construyo a continuación:

library(randomForest)

raw <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/credit-screening/crx.data",
    sep = ",", na.strings = "?")

dat <- raw
dat$V14 <- dat$V6 <- NULL    # me da igual
dat <- na.omit(dat)          # ídem

modelo <- randomForest(V16 ~ ., data = dat)

Fijémonos en el sujeto 100, a quien se le deniega el crédito (suponiendo, ¡mal hecho!, que el punto de corte de la probabilidad para concederlo es el 50%), y la variable $V8$. Podemos ver cuál sería el score del cliente modificando esa variable entre su valor real y el máximo del rango dejando las demás tal cual:

El telón de Aquiles del big data es el sesgo. Me gustaría hablar más de ello, pero me agarra de la pluma uno de esos NDAs. Así que hablaré de otra cosa.

Si le preguntas a la gente cuántos hermanos son en la familia, el promedio del resultado tenderá a ser superior al número medio de hijos por familia. Esencialmente, porque no estás muestreando familias sino hijos. El tautológico hecho de que las familias con más hijos tengan más hijos hace que estén sobrerrepresentadas en la muestra.

Es

$$ \sum_i \Phi(y_i, f_1(x_i)) > \sum_i \Phi(y_i, f_1(x_i) - \lambda \nabla \Phi(y_i, f_1(x_i)) \sim$$ $$ \sim \sum_i \Phi(y_i, f_1(x_i) - \lambda f_2(x_i))$$

Por supuesto, el lector se preguntará muchas cosas, entre las que destaco:

• ¿Qué representa cada uno de los elementos que aparecen en la línea anterior?
• ¿Qué parte de ella es solo casi siempre cierta?
• ¿Qué tiene todo eso que ver con GBM?

… (por motivos que importan pero no debo revelar a mis lectores) aprovecho para criticar a esos tipos que, vistiendo como yo, insisten reiteradamente a sus analistas en que les proporcionen un número. Un número que tiene que ser cerrado, indiscutible, pivotal.

A esos que gastan traje y corbata como yo hoy les horroriza la varianza. Le espantan, seguro, esos punticos que tan opotunamente coloca Kiko Llaneras alrededor de las medias de este estupendo

Que es una manera de matizar

Porque, recordemos,

• no solo con datos tomamos decisiones informadas: las prioris (experiencia cuantificada) tienen su importancia
• no podemos obtener datos que justifiquen todas, todas, todas las decisiones.

y <- c(rnorm(1000), rnorm(2000, 1, 0.5))

es una mezcla de dos normales (N(0, 1) y N(1, 0.5)) con pesos 1/3 y 2/3 respectivamente. Pero, ¿cómo podríamos estimar los parámetros a partir de esos datos?

Se puede usar, p.e., flexmix, que implementa eso del EM. Pero en el librillo de este maestrillo dice

library(rstan)

y <- c(rnorm(1000), rnorm(2000, 1, 0.5))

codigo <- "
data {
  int<lower=1> K; // number of mixture components
  int<lower=1> N; // number of data points
  real y[N]; // observations
}

parameters {
  simplex[K] theta; // mixing proportions
  real mu[K]; // locations of mixture components
  real<lower=0> sigma[K]; // scales of mixture components
}

model {
  real ps[K]; // temp for log component densities

  sigma ~ cauchy(0,2.5);
  mu ~ normal(0,10);

  for (n in 1:N) {
    for (k in 1:K) {
      ps[k] <- log(theta[k]) + normal_log(y[n],mu[k], sigma[k]);
    }
    increment_log_prob(log_sum_exp(ps));
  }
}"

fit <- stan(model_code = codigo,
            data = list(K = 2, N = length(y), y = y),
            iter=48000, warmup=2000,
            chains=1, thin=10)

En el código anterior no sé si queda claro cómo cada punto $latex y_i$ sigue una distribución (condicionada a los parámetros) con densidad $latex \theta_1 \phi(y_i, \mu_1, \sigma_1) + \theta_2 \phi(y_i, \mu_2, \sigma_2)$.