Estadística Bayesiana

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.

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í:

Porque ahí puedes tomarte una foto tal que

o

y luego tuitear cosas como

After #StrataHadoop: @jayusor and me in front of Bayes grave with my new @OReillyMedia book (emulating @gilbellosta) pic.twitter.com/SlguJIeLw0

— Antonio Sánchez Chinchón (@aschinchon) June 4, 2016

Para mayor referencia (y por tenerlo a mano cuando vuelva),

Las diapositivas de mi charla Datos, modelos y parámetros en el grupo Machine Learning Spain pueden verse/bajarse de aquí.

Porque voy a dar una charla en él. Es este jueves, por la tarde, en el Campus de Google de Madrid (los detalles).

Se tratará de una introducción a y justificación de aproximaciones más bayesianas de lo habitual a problemas reales del análisis de datos. Que comenzará con una explicación sobre cuándo 100% no significa 100% para terminar con lo que viene siéndome habitual últimamente: un ejemplo en rstan con su discusión.

Hago pública por su interés (parte de) una respuesta de Ramón Díaz Uriarte a un correo mío en el que yo sugería

que una vez que sabes especificar un modelo probabilístico para unos datos, p.e.,

• para la regresión lineal, y ~ N(a0 + a1 x1 +..., sigma)),
• para el test de Student, y0 ~ N(mu, sigma); y1 ~ N(mu + delta, sigma),
• etc. no hace falta saber qué es lm, ni el test de Student, ni nada. Cero teoría; sobre todo, de teoría tipo recetario. Se especifica el modelo (con una determinada sintaxis), se deja correr la cosa y a interpretar.

Su respuesta:

Este jueves (2016-02-11), a las 19:00, hablaré de rstan y de rstanarm en Medialab-Prado dentro de la reunión de usuarios de R de Madrid. Con el concurso de estos paquetes, replantearé tres problemas estadísticos conocidos desde una óptica bayesiana:

• Pruebas de hipótesis
• Regresión lineal
• Modelos estructurales de series temporales

Si quieres asistir, reserva tu plaza aquí.

Probablemente, discutiré todos esos modelos en estas páginas en los próximos días, además de colgar las diapositivas y sus fuentes.

Y termino con lo de los intervalos. Me refiero a esto y esto.

Nunca me habría atrevido a escribir sobre el tema, y exponerme, de paso, a la muy razonadas explicaciones de quienes tuvieron a bien comentarlas, si no hubiese sido por un tema personal: el recuerdo de la frustración que me supuso hacerme en su día con la teoría subyacente tanto a las pruebas de hipótesis como a la construcción de intervalos de confianza.

A veces se hacen encuestas sobre temas sobre los que los encuestados son reticentes a revelar la verdad (p.e., ¿es Vd. un zombi?). Un procedimiento conocido para recabar tal tipo de información es el siguiente:

• Se le invita al encuestado a tirar al aire una moneda con las caras etiquetadas con sí y no; la moneda no es una moneda porque tiene una probabidad conocida (y distinta del 50%) de caer en sí.
• El encuestado responde sí si la respuesta a la pregunta y el resultado de la tirada de la moneda coinciden y no en caso contrario.

A partir de la proporción de respuestas positivas y conocida la probabilidad del sí de la moneda, $latex q$, es posible estimar la proporción $latex \theta$ de respuestas positivas a la pregunta de subyacente de interés en la muestra. Efectivamente, los síes tienen una distribución binomial $latex B(p) = B(q\theta + (1-q)(1-\theta))$ y, una vez estimado (por máxima verosimilitud) $latex \hat{p}$, puede despejarse $latex \hat{p}$ de $latex \hat{p} = q\hat{\theta} + (1-q)(1-\hat{\theta})$ para obtener

El problema en cuestión, que se ve, surgió durante la II Guerra Mundial, es el siguiente: se capturan tanques del enemigo y se anotan los números de serie, supuestos sucesivos. ¿Cuál es la mejor estimación del número total de tanques fabricados por el enemigo?

Si se capturan k, la distribución del máximo número observado, m, en función del número no observado (nuestro parámetro) de tanques es

$$ f(N;m,k)=\frac{\binom{m-1}{k-1}}{\binom{N}{k}}$$

y como esta función es decreciente en $latex N$, la estimación por máxima verosimilitud es $latex \hat{N} = m$.