Estou tentando calcular o intervalo de 95% credível da seguinte distribuição posterior. Não consegui encontrar a função em R para isso, mas a abordagem abaixo está correta?
x <- seq(0.4,12,0.4)
px <- c(0,0, 0, 0, 0, 0, 0.0002, 0.0037, 0.018, 0.06, 0.22 ,0.43, 0.64,0.7579, 0.7870, 0.72, 0.555, 0.37, 0.24, 0.11, 0.07, 0.02, 0.009, 0.005, 0.0001, 0,0.0002, 0, 0, 0)
plot(x,px, type="l")
mm <- sum(x*px)/sum(px)
var <- (sum((x)^2*px)/sum(px)) - (mm^2)
cat("95% credible interval: ", round(mm -1.96*sqrt(var),3), "-", round(mm + 1.96*sqrt(var),3),"\n")
bayesian
descriptive-statistics
credible-interval
user19758
fonte
fonte
Respostas:
Conforme observado por Henry , você está assumindo a distribuição normal e não há problema se seus dados seguirem a distribuição normal, mas estarão incorretos se você não puder assumir a distribuição normal para eles. Abaixo, descrevo duas abordagens diferentes que você pode usar para distribuição desconhecida, considerando apenas pontos de dados
x
e estimativas de densidade correspondentespx
.A primeira coisa a considerar é o que exatamente você deseja resumir usando seus intervalos. Por exemplo, você pode estar interessado nos intervalos obtidos usando quantis, mas também na região de maior densidade (veja aqui ou aqui ) da sua distribuição. Embora isso não deva fazer muita diferença (se houver) em casos simples, como distribuições simétricas e unimodais, isso fará diferença para distribuições mais "complicadas". Geralmente, os quantis fornecem um intervalo contendo massa de probabilidade concentrada em torno da mediana (os médios da sua distribuição), enquanto a região de maior densidade é uma região em torno dos modos100 α% da distribuição. Isso ficará mais claro se você comparar as duas parcelas da figura abaixo - os quantis "cortam" a distribuição verticalmente, enquanto a região de maior densidade "corta" horizontalmente.
A próxima coisa a considerar é como lidar com o fato de você ter informações incompletas sobre a distribuição (assumindo que estamos falando de distribuição contínua, você tem apenas alguns pontos e não uma função). O que você pode fazer é pegar os valores "como estão" ou usar algum tipo de interpolação ou suavização para obter os valores "entre".
Uma abordagem seria usar a interpolação linear (veja
?approxfun
em R) ou, alternativamente, algo mais suave como splines (veja?splinefun
em R). Se você escolher essa abordagem, lembre-se de que os algoritmos de interpolação não têm conhecimento de domínio sobre seus dados e podem retornar resultados inválidos, como valores abaixo de zero, etc.A segunda abordagem que você pode considerar é usar a distribuição de densidade / mistura do kernel para aproximar sua distribuição usando os dados que você possui. A parte complicada aqui é decidir sobre a largura de banda ideal.
Em seguida, você encontrará os intervalos de interesse. Você pode prosseguir numericamente ou por simulação.
1a) Amostragem para obter intervalos quantílicos
1b) Amostragem para obter a região de maior densidade
2a) Encontre quantis numericamente
2b) Encontre a região de maior densidade numericamente
Como você pode ver nas plotagens abaixo, no caso de distribuição simétrica unimodal, ambos os métodos retornam o mesmo intervalo.
fonte