Era como
y se ha convertido en
¡Qué horror!
Coda: En otra página de la Wikipedia en la que he caído después por azar he leído la siguiente frase (que por algún motivo encuentro relevante insertar aquí):
Los ríos arrastran sedimentos que consiguen colmatar y rellenar de lodo los lagos. Además, la proliferación de ciertas plantas, como el lirio acuático, los obstruye por completo.
Más sobre lo de ayer. O más bien, una justificación por analogía.
Con monedas.
Tiras una moneda 100 veces y obtienes 60 caras. Tienes una priori $latex B(a,b)$ (beta). Tomas una muestra de valores $latex p_i$ con esa distribución y para cada una de ellas repites el experimento, es decir, obtienes lo que en R se expresaría de la forma
rbinom(1, 100, p[i])Si te quedas los valores $p_i$ tales que esa simulación es 60, enhorabuena, tienes una muestra de la distribución a posteriori.
Tengo el sistema
$$ m = \frac{a}{a+b}$$ $$ v = \frac{ab}{(a+b)^2 (a+b+1)}$$
en los que alguien descubrirá cosas relativas a la distribución beta.
Interesa despejar $latex a$ y $latex b$. Pero solo soy un exmatemático perezoso, disléxico y con déficit de tiempo y atención. Así que tacacata y…
$$ a = \frac{-m^3 + m^2 - mv}{v}$$ $$ b = \frac{m^3 - 2m^2 + mv + m -v}{v}$$
Aquí, contracorriente. Dejamos aparcado el big data y le damos a lo que nos da de comer. Entre otras cosas, este pequeño experimento con muy pequeños datos (¿tres?).
La aplicación es real. Y los datos pequeños porque son carísimos.
Se puede suponer que tienen distribución beta de parámetros desconocidos. Nos interesa la media muestral de unas pocas observaciones: dos, tres, cuatro,… En particular, qué distribución tiene.
Si fuesen muchos, podríamos aplicar el teorema central del límite (que funciona estupendamente incluso con valores no muy grandes). Pero la suma de pocas observaciones beta no tiene una distribución con nombre (que yo sepa). Pero podemos usar un viejo truco (parecido al de la aproximación de Welch para el número de grados de libertad de la prueba de Student cuando las varianzas son desiguales):
La estadística bayesiana se enseña en cursos de estadística (y, frecuentemente, envuelto en un aparataje matemático tan ofuscante como innecesario). Lo malo es que en los cursos y textos de estadística no existe información previa. La información previa sobre los fenómenos en los que se utilizaría la estadística bayesiana están en las aplicaciones, extramuros del muy agnóstico mundo de la estadística y la matemática.
Por eso, a los autores de los libros de estadística bayesiana y quienes enseñan cursos sobre lo mismo, enfrentados al problema de llenar de sentido la problemática distribución a priori, no se les ocurre nada mejor que discutir muy sesudamente la excepción (la priori no informativa) en lugar de la regla (la priori informativa). Reto al lector escéptico a que repase cualquier manual en la materia (que no haya sido escrito por Gelman) y compare el espacio que dedican a la selección de prioris no informativas con el de convenir una priori informativa decente.
A partir de los comentarios de Olivier Núñez a mi entrada anterior casi homónima, se nos ha ocurrido a ambos de forma independiente y simultánea una manera alternativa de calcular el intervalo: minimizando su longitud.
a <- 3
b <- 5
alfa <- 0.05
# versión de la entrada anterior:
f <- function(x){
(dbeta(x[2], a, b) - dbeta(x[1], a, b))^2 +
(pbeta(x[2], a, b) - pbeta(x[1], a, b) -1 + alfa)^2
}
res <- optim(c(a/(a+b), a/(a+b)), f)
res$par
#[1] 0.08052535 0.68463436
# nueva versión
f.alt <- function(x){
qbeta(x+0.95, a, b) - qbeta(x, a, b)
}
res.alt <- optim(0.025, f.alt)
qbeta(c(res.alt$par, res.alt$par + 0.95), a, b)
#[1] 0.08054388 0.68464900Tengo un parámetro, la p de una binomial, que supongo distribuido según una beta. Me da igual para el caso si la distribución a priori es o no informativa. Solo digo que la distribución a posteriori es otra beta con parámetros a y b.
Quiero construir un intervalo de credibilidad para p, es decir, encontrar un subintervalo de [0,1]
Gráficamente,
