Tengo 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]
- dentro del cual la densidad de la beta sea mayor que fuera y que
- capture
de la probabilidad total.
Gráficamente,
Y he aquí el código:
a <- 3 b <- 5 alfa <- 0.05 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) x <- 1:100 / 100 plot(x, dbeta(x, a, b), type = "l", ylab = "densidad") lines(c(res$par[1], res$par[1]), c(0, dbeta(res$par[1], a, b)), col = "red") lines(c(res$par[2], res$par[2]), c(0, dbeta(res$par[2], a, b)), col = "red") lines(c(res$par[1], res$par[2]), rep(dbeta(res$par[2], a, b), 2), col = "red")
La función que se optimiza tiene como argumentos los puntos inicial y final del intervalo y penaliza:
- Que la densidad en dichos punto sea distinta.
- Que la suma de las probabilidades de las colas descartadas sea distinta de
.
Una posible mejora en el código anterior sería pasarle a optim
mejores puntos de partida: en lugar de la media de la distribución para ambos casos, la media más (y menos) dos veces la desviación estándar.
Porqué no utilizar los cuantiles de la beta?
lim_inf = qbeta(alfa/2,a,b) ; lim_sup = qbeta(1-alfa/2,a,b)
Ya entiendo: tu método minimiza la longitud del intervalo. Pues entonces, estos cuantiles pueden servir al menos de puntos de partida.