A menudo se da el caso de que un intervalo de confianza con una cobertura del 95% es muy similar a un intervalo creíble que contiene el 95% de la densidad posterior. Esto sucede cuando el prior es uniforme o casi uniforme en el último caso. Por lo tanto, un intervalo de confianza a menudo se puede utilizar para aproximar un intervalo creíble y viceversa. Es importante destacar que podemos concluir de esto que la interpretación errónea muy difamada de un intervalo de confianza como un intervalo creíble tiene poca o ninguna importancia práctica para muchos casos de uso simple.
Hay una serie de ejemplos de casos en los que esto no sucede, sin embargo, todos parecen ser engañados por los defensores de las estadísticas bayesianas en un intento de demostrar que hay algo mal con el enfoque frecuentista. En estos ejemplos, vemos que el intervalo de confianza contiene valores imposibles, etc., que se supone que muestran que no tienen sentido.
No quiero volver sobre esos ejemplos, o una discusión filosófica de Bayesiano vs Frequentista.
Solo estoy buscando ejemplos de lo contrario. ¿Existen casos en los que la confianza y los intervalos creíbles son sustancialmente diferentes, y el intervalo proporcionado por el procedimiento de confianza es claramente superior?
Para aclarar: se trata de la situación en la que generalmente se espera que el intervalo creíble coincida con el intervalo de confianza correspondiente, es decir, cuando se usan anteriores planos, uniformes, etc. No estoy interesado en el caso en el que alguien elige un mal arbitrario anterior.
EDITAR: En respuesta a la respuesta de @JaeHyeok Shin a continuación, debo estar en desacuerdo con que su ejemplo utiliza la probabilidad correcta. Utilicé el cálculo bayesiano aproximado para estimar la distribución posterior correcta para theta a continuación en R:
### Methods ###
# Packages
require(HDInterval)
# Define the likelihood
like <- function(k = 1.2, theta = 0, n_print = 1e5){
x = NULL
rule = FALSE
while(!rule){
x = c(x, rnorm(1, theta, 1))
n = length(x)
x_bar = mean(x)
rule = sqrt(n)*abs(x_bar) > k
if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
}
return(x)
}
# Plot results
plot_res <- function(chain, i){
par(mfrow = c(2, 1))
plot(chain[1:i, 1], type = "l", ylab = "Theta", panel.first = grid())
hist(chain[1:i, 1], breaks = 20, col = "Grey", main = "", xlab = "Theta")
}
### Generate target data ###
set.seed(0123)
X = like(theta = 0)
m = mean(X)
### Get posterior estimate of theta via ABC ###
tol = list(m = 1)
nBurn = 1e3
nStep = 1e4
# Initialize MCMC chain
chain = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = c("theta", "mean")
chain$theta[1] = rnorm(1, 0, 10)
# Run ABC
for(i in 2:nStep){
theta = rnorm(1, chain[i - 1, 1], 10)
prop = like(theta = theta)
m_prop = mean(prop)
if(abs(m_prop - m) < tol$m){
chain[i,] = c(theta, m_prop)
}else{
chain[i, ] = chain[i - 1, ]
}
if(i %% 100 == 0){
print(paste0(i, "/", nStep))
plot_res(chain, i)
}
}
# Remove burn-in
chain = chain[-(1:nBurn), ]
# Results
plot_res(chain, nrow(chain))
as.numeric(hdi(chain[, 1], credMass = 0.95))
Este es el intervalo de 95% creíble:
> as.numeric(hdi(chain[, 1], credMass = 0.95))
[1] -1.400304 1.527371
EDITAR # 2:
Aquí hay una actualización después de los comentarios de @JaeHyeok Shin. Estoy tratando de mantenerlo lo más simple posible, pero el guión se volvió un poco más complicado. Cambios principales:
- Ahora usando una tolerancia de 0.001 para la media (era 1)
- Mayor número de pasos a 500k para tener en cuenta una tolerancia menor
- Disminuyó el SD de la distribución de la propuesta a 1 para tener en cuenta una tolerancia menor (era 10)
- Se agregó la probabilidad de tormenta simple con n = 2k para la comparación
- Se agregó el tamaño de muestra (n) como estadística de resumen, establezca la tolerancia en 0.5 * n_target
Aquí está el código:
### Methods ###
# Packages
require(HDInterval)
# Define the likelihood
like <- function(k = 1.3, theta = 0, n_print = 1e5, n_max = Inf){
x = NULL
rule = FALSE
while(!rule){
x = c(x, rnorm(1, theta, 1))
n = length(x)
x_bar = mean(x)
rule = sqrt(n)*abs(x_bar) > k
if(!rule){
rule = ifelse(n > n_max, TRUE, FALSE)
}
if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
}
return(x)
}
# Define the likelihood 2
like2 <- function(theta = 0, n){
x = rnorm(n, theta, 1)
return(x)
}
# Plot results
plot_res <- function(chain, chain2, i, main = ""){
par(mfrow = c(2, 2))
plot(chain[1:i, 1], type = "l", ylab = "Theta", main = "Chain 1", panel.first = grid())
hist(chain[1:i, 1], breaks = 20, col = "Grey", main = main, xlab = "Theta")
plot(chain2[1:i, 1], type = "l", ylab = "Theta", main = "Chain 2", panel.first = grid())
hist(chain2[1:i, 1], breaks = 20, col = "Grey", main = main, xlab = "Theta")
}
### Generate target data ###
set.seed(01234)
X = like(theta = 0, n_print = 1e5, n_max = 1e15)
m = mean(X)
n = length(X)
main = c(paste0("target mean = ", round(m, 3)), paste0("target n = ", n))
### Get posterior estimate of theta via ABC ###
tol = list(m = .001, n = .5*n)
nBurn = 1e3
nStep = 5e5
# Initialize MCMC chain
chain = chain2 = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = colnames(chain2) = c("theta", "mean")
chain$theta[1] = chain2$theta[1] = rnorm(1, 0, 1)
# Run ABC
for(i in 2:nStep){
# Chain 1
theta1 = rnorm(1, chain[i - 1, 1], 1)
prop = like(theta = theta1, n_max = n*(1 + tol$n))
m_prop = mean(prop)
n_prop = length(prop)
if(abs(m_prop - m) < tol$m &&
abs(n_prop - n) < tol$n){
chain[i,] = c(theta1, m_prop)
}else{
chain[i, ] = chain[i - 1, ]
}
# Chain 2
theta2 = rnorm(1, chain2[i - 1, 1], 1)
prop2 = like2(theta = theta2, n = 2000)
m_prop2 = mean(prop2)
if(abs(m_prop2 - m) < tol$m){
chain2[i,] = c(theta2, m_prop2)
}else{
chain2[i, ] = chain2[i - 1, ]
}
if(i %% 1e3 == 0){
print(paste0(i, "/", nStep))
plot_res(chain, chain2, i, main = main)
}
}
# Remove burn-in
nBurn = max(which(is.na(chain$mean) | is.na(chain2$mean)))
chain = chain[ -(1:nBurn), ]
chain2 = chain2[-(1:nBurn), ]
# Results
plot_res(chain, chain2, nrow(chain), main = main)
hdi1 = as.numeric(hdi(chain[, 1], credMass = 0.95))
hdi2 = as.numeric(hdi(chain2[, 1], credMass = 0.95))
2*1.96/sqrt(2e3)
diff(hdi1)
diff(hdi2)
Los resultados, donde hdi1 es mi "probabilidad" y hdi2 es la simple rnorm (n, theta, 1):
> 2*1.96/sqrt(2e3)
[1] 0.08765386
> diff(hdi1)
[1] 1.087125
> diff(hdi2)
[1] 0.07499163
Entonces, después de reducir la tolerancia lo suficiente, y a expensas de muchos más pasos de MCMC, podemos ver el ancho CrI esperado para el modelo rnorm.
Respuestas:
En un diseño experimental secuencial, el intervalo creíble puede ser engañoso.
(Descargo de responsabilidad: no estoy argumentando que no es razonable, es perfectamente razonable en el razonamiento bayesiano y no es engañoso desde la perspectiva del punto de vista bayesiano).
Mensaje para llevar a casa: si está interesado en tener una garantía frecuentista, debe tener cuidado al usar las herramientas de inferencia bayesianas, que siempre son válidas para las garantías bayesianas, pero no siempre para las frecuentes.
(Aprendí este ejemplo de la increíble conferencia de Larry. Esta nota contiene muchas discusiones interesantes sobre la sutil diferencia entre los marcos frecuentista y bayesiano. Http://www.stat.cmu.edu/~larry/=stat705/Lecture14.pdf )
EDITAR En el ABC de Livid, el valor de tolerancia es demasiado grande, por lo que incluso para la configuración estándar donde muestreamos un número fijo de observaciones, no proporciona un CR correcto. No estoy familiarizado con ABC, pero si solo cambié el valor de tol a 0.05, podemos tener un CR muy sesgado de la siguiente manera
> as.numeric(hdi(chain[, 1], credMass = 0.95)) [1] -0.01711265 0.14253673
Por supuesto, la cadena no está bien estabilizada, pero incluso si aumentamos la longitud de la cadena, podemos obtener un CR similar, sesgado a una parte positiva.
fuente
Dado que el intervalo creíble se forma a partir de la distribución posterior, en base a una distribución previa estipulada, puede construir fácilmente un intervalo creíble muy malo utilizando una distribución previa que esté muy concentrada en valores de parámetros altamente inverosímiles. Puede hacer un intervalo creíble que no tenga "sentido" mediante el uso de una distribución previa que esté completamente concentrada en valores de parámetros imposibles .
fuente
Si estamos usando un previo plano, este es simplemente un juego en el que tratamos de llegar a un previo plano en una reparametrización que no tiene sentido.
Esta es la razón por la cual muchos bayesianos se oponen a las anteriores planas.
fuente