Geralmente, um intervalo de confiança com cobertura de 95% é muito semelhante a um intervalo credível que contém 95% da densidade posterior. Isso acontece quando o prior é uniforme ou quase uniforme no último caso. Assim, um intervalo de confiança geralmente pode ser usado para aproximar um intervalo confiável e vice-versa. É importante ressaltar que podemos concluir que a má interpretação muito difamada de um intervalo de confiança como um intervalo credível tem pouca ou nenhuma importância prática para muitos casos de uso simples.
Existem vários exemplos de casos em que isso não acontece, no entanto, todos parecem ser escolhidos pelos defensores das estatísticas bayesianas, na tentativa de provar que há algo errado com a abordagem freqüentista. Nestes exemplos, vemos que o intervalo de confiança contém valores impossíveis, etc., o que deve mostrar que eles não fazem sentido.
Não quero voltar atrás nesses exemplos ou em uma discussão filosófica sobre Bayesiana vs Frequentista.
Estou apenas procurando exemplos do oposto. Existem casos em que os intervalos de confiança e credibilidade são substancialmente diferentes e o intervalo fornecido pelo procedimento de confiança é claramente superior?
Para esclarecer: Trata-se da situação em que geralmente se espera que o intervalo credível coincida com o intervalo de confiança correspondente, ou seja, ao usar anteriores simples, uniformes, etc. Não estou interessado no caso em que alguém escolhe um prioritário arbitrariamente ruim.
EDIT: Em resposta à resposta de @JaeHyeok Shin abaixo, devo discordar de que o exemplo dele usa a probabilidade correta. Eu usei o cálculo bayesiano aproximado para estimar a distribuição posterior correta para teta abaixo em 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 é o intervalo credível de 95%:
> as.numeric(hdi(chain[, 1], credMass = 0.95))
[1] -1.400304 1.527371
EDIT # 2:
Aqui está uma atualização após os comentários de @JaeHyeok Shin. Estou tentando mantê-lo o mais simples possível, mas o script ficou um pouco mais complicado. Principais mudanças:
- Agora, usando uma tolerância de 0,001 para a média (era 1)
- Aumento do número de etapas para 500k, para diminuir a tolerância
- Diminuiu o SD da distribuição da proposta para 1, para levar em conta a menor tolerância (era 10)
- Adicionada a probabilidade rnorm simples com n = 2k para comparação
- Adicionado o tamanho da amostra (n) como uma estatística resumida, defina a tolerância como 0,5 * n_target
Aqui está o 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)
Os resultados, em que hdi1 é minha "probabilidade" e hdi2 é o rnorm simples (n, theta, 1):
> 2*1.96/sqrt(2e3)
[1] 0.08765386
> diff(hdi1)
[1] 1.087125
> diff(hdi2)
[1] 0.07499163
Portanto, depois de diminuir a tolerância suficientemente e às custas de muitas outras etapas do MCMC, podemos ver a largura de CrI esperada para o modelo rnorm.
Respostas:
Em um projeto experimental seqüencial, o intervalo credível pode ser enganoso.
(Isenção de responsabilidade: não estou argumentando que não é razoável - é perfeitamente razoável no raciocínio bayesiano e não é enganoso na perspectiva do ponto de vista bayesiano.)
Para um exemplo simples, digamos que temos uma máquina que nos fornece uma amostra aleatóriaX de N(θ,1) com θ desconhecido . Em vez de desenhar n iid amostras, extraímos amostras até n−−√X¯n>k k N N=inf{n≥1:n−−√X¯n>k}.
Leve a mensagem para casa: Se você estiver interessado em ter uma garantia freqüentista, tenha cuidado ao usar ferramentas de inferência bayesiana que são sempre válidas para garantias bayesianas, mas nem sempre para garantias freqüentadoras.
(Aprendi esse exemplo com a impressionante palestra de Larry. Esta nota contém muitas discussões interessantes sobre a diferença sutil entre estruturas freqüentistas e bayesianas. Http://www.stat.cmu.edu/~larry/=stat705/Lecture14.pdf )
EDIT No ABC da Livid, o valor da tolerância é muito grande, portanto, mesmo para a configuração padrão em que amostramos um número fixo de observações, ele não fornece um CR correto. Eu não estou familiarizado com o ABC, mas se eu apenas alterei o valor de tol para 0,05, podemos ter um CR muito distorcido da seguinte maneira
> as.numeric(hdi(chain[, 1], credMass = 0.95)) [1] -0.01711265 0.14253673
Obviamente, a corrente não está bem estabilizada, mas mesmo se aumentarmos o comprimento da corrente, podemos obter um CR semelhante - inclinado para uma parte positiva.
fonte
Como o intervalo credível é formado a partir da distribuição posterior, com base em uma distribuição anterior estipulada, você pode facilmente construir um intervalo credível muito ruim usando uma distribuição anterior muito concentrada em valores de parâmetros altamente implausíveis. Você pode criar um intervalo credível que não "faça sentido" usando uma distribuição anterior totalmente concentrada em valores de parâmetros impossíveis .
fonte
Se estamos usando um flat anterior, este é simplesmente um jogo em que tentamos criar um flat anterior em uma reparameterização que não faz sentido.
É por isso que muitos bayesianos se opõem a priores planos.
fonte