Estadística

El otro día asistí a la enésima confusión sobre k-medias y k-vecinos. Que lo es, más en general, sobre el clústering contra modelos locales de la clase que sean, desde k-vecinos hasta el filtrado colaborativo. Veamos si esta comparación que traigo hoy a mis páginas contribuye a erradicar dicha confusión.

k-medias es como las elecciones. Hace poco tuvimos unas en España. Alguien decidió (aproximadamente) que k = 4 y nos pidió, a nosotros, punticos del espacio, identificar el centroide más próximo a nosotros para que lo votásemos. Pues eso, la misma frustración que muchos dizque sintieron teniendo que elegir entre partidos/centroides subjetivamente igual de alejados de los intereses de uno es la que sienten nuestros punticos cuando los procrusteamos para asociarlos al totum revolutum de los clientes estrella, etc.

Esto que cuento hoy puede ser muy útil: cómo mejorar los GLMs mediante la introducción de prioris (casi) a voluntad sobre los coeficientes. Usando el paquete arm de R, claro.

De momento y porque aún tengo sucios los datos sobre los que me gustaría aplicar el modelo, extraeré un ejemplo de la ayuda de la función principal del paquete, bayesglm.

Primero, preparo unos datos:

n <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1
b1 <- 1.5
b2 <- 2
y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2))

Comenzamos con un glm de toda la vida.

Continuando con la entrada anterior, ahora, números.

Primero, el planteamiento (cuatro partidos, etc.):

probs <- c(4, 3, 2, 1)
probs <- probs / sum(probs)
partidos <- letters[1:length(probs)]

Nos hará falta más adelante

library(plyr)
library(rstan)
library(ggplot2)
library(reshape2)

Sigo con el proceso de muestreo. Reitero: cada encuestador enseña al encuestado una tarjeta al azar donde aparece el nombre de dos partidos y le pregunta si ha votado (o piensa votar) a alguno de ellos.

n <- 3000
resultados <- data.frame(
  tarjeta = sample(1:nrow(tarjetas), n, replace = T),
  partido = sample(partidos, n, prob = probs, replace = T))
resultados <- data.frame(
  tarjetas[resultados$tarjeta,],
  partido = resultados$partido)
resultados$coincide <- resultados$partido == resultados$partido1 |
  resultados$partido == resultados$partido2

# proporciones reales en la muestra
props.muestra <- table(resultados$partido) / nrow(resultados)

# resultados agregados (por tarjeta)
resultados.agg <- ddply(
    resultados, .(partido1, partido2),
    summarize,
    total = length(partido1),
    coincidencias = sum(coincide))

Y

Lo de que la gente que miente al ser encuestada se ha esgrimido frecuentemente en los últimos días. Inspirado en esto, se me ha ocurrido (posiblemente reocurrido: es fácil que la idea sea conocida, sobre todo si resulta ser buena) el siguiente procedimiento para la realización de encuestas electorales.

• El encuestador va provisto de una colección de cartulinas en las que aparecen parejas de nombres de partidos políticos.
• El encuestador muestra al encuestado una cartulina al azar dentro de su colección.
• El encuestador pregunta al encuestado si ha votado (o piensa votar) a alguno de ellos.
• Se registran los partidos mostrados y la respuesta, positiva o negativa, del encuestado.

Con una versión del procedimiento que describo en la entrada que enlazo más arriba, se podrían redescubrir las opciones de la población subyacente, aun ignorando el de cada uno de los encuestados. No sé cuál sería (si no se me adelanta nadie, igual la hago yo) el procedimiento, pero seguro que no es tan complicado como para que Stan no pueda con ello.

Me piden que opine sobre lo de las encuestas electorales y su error. Vaya por delante mi confesión de que de eso sé poco. Soy matemático, no estadístico, y uno de los mayores huecos (¿simas?) de mi formación estadística tiene que ver con todo lo relativo al muestreo. Así que, con la valentía que aporta la ignorancia, procedo.

El primer gran problema con las encuestas electorales es que confunden países con urnas y gente con bolas de colores. Si en una urna hay N bolas (de colores distintos) y queremos estimar su número mediante una extracción de n bolas, existe un margen de error debido a que en lugar de ver todos los datos uno ve únicamente una muestra.

Liberados del estrecho ámbito de nuestra original mentira sugerente gracias a la relación que descubrimos entre residuos y gradientes cuando las pérdidas son cuadráticas podemos adentrarnos en ámbitos más extensos.

Lo que discutimos del gradiente tiene una interpretación fácilmente inteligible en el caso de pérdidas cuadráticas. Pero ni la pérdida de interpretabilidad nos impide extender el razonamiento de la entrada anterior a funciones de pérdida distintas de la cuadrática siempre que podamos calcular un gradiente.

Para minimizar una función $latex \phi(x)$ es habitual utilizar un procedimiento iterativo: a partir de un punto inicial $latex x_0$ se salta a $latex x_1 = x_0 - \lambda \nabla \phi(x_0)$ (donde $latex \lambda$ es un número pequeño predefinido de antemano), y de ahí, sucesivamente, a

$$ x_n = x_{n-1} - \lambda \nabla \phi(x_{n-1}).$$

Porque, típicamente, como cuando uno está en el monte y da un paso corto en la dirección opuesta a la de máxima pendiente, sucede que

Hace un tiempo resumí los GBMs (Gradient Boosting Machines) en una línea. Hoy comienzo una serie de varias entradas para que nadie tenga excusa de no saber de qué va la cosa. Arranco con una mentira sugerente. Porque lo que voy a contar no es del todo cierto, pero motiva lo que vendrá después.

Consideremos un conjunto de datos medio famoso: el de los precios de los alquileres en Múchich. Comencemos con un modelo sencillo, una regresión lineal que relacione el precio del alquiler con los metros cuadrados, i.e.,

Tengo la sensación de que un lenguaje funcional (como Scala) está particularmente bien adaptado al tipo de operaciones que exige MCMC.

Juzguen Vds.

Primero, genero datos en R:

datos <- rnorm(500, 0.7, 1)
writeLines(as.character(datos), "/tmp/datos.txt")

Son de una normal con media 0.7. En el modelo que vamos a crear, suponemos conocida (e igual a 1) la varianza de la normal y trataremos de estimar la media suponiéndole una distribución a priori normal estándar. Y con Scala, así:

Aunque esta entrada es sin duda resabida de los más de mis lectores, quedarán los que aún no sepan que ciertas distribuciones no tienen media. Condición necesaria para que una distribución la tenga es que

$$ \int_{-\infty}^\infty |x| f(x) dx$$

tenga un valor finito, cosa que, por ejemplo, no cumple la de Cauchy. Igual hay a quien esto le parece una rareza matemática, un entretenimiento de math kiddies sin implicaciones prácticas. Además, porque para que que la integral anterior diverja se necesita que las distribuciones puedan tomar valores arbitrariamente altos y las que se manejan en la práctica están acotadas si no por el número de átomos del universo por el de céntimos de bolívar venezolano necesarios para comprar todas las cosas que caben en el ancho mundo.