Cómo promediar varias distribuciones posteriors a partir de una simulación Monte Carlo

Question

Digamos que usted produce varias distribuciones posteriors a partir de diferentes ejecuciones del mismo modelo bajo diferentes semillas. Es decir, usted tiene algo como lo siguiente:

library(ggplot2)
set.seed(10100101)
N <- 150;R <- 25
normals <- lapply(1:R,function(x){rnorm(n = N, mean = 0, sd   = 1)})
df <- data.frame(Reduce(cbind, normals))
colnames(df) <- 1:ncol(df)
df2 <- stack(df)
ggplot(df2, aes(values, fill = ind)) + geom_density(alpha = 0.000001) + theme(legend.position="none")

Me gustaría añadir al gráfico la "media de las posteriores" producida por las diferentes ejecuciones de mi modelo. Al principio pensé que sería tan fácil como aplicar la media de las posteriores. Sin embargo, estoy obteniendo algo que converge a cero porque simplemente estoy promediando valores que provienen de una distribución normal estándar como se puede ver en la imagen de abajo:

library(ggplot2)
set.seed(10100101)
N <- 150;R <- 25
normals <- lapply(1:R,function(x){rnorm(n = N, mean = 0, sd   = 1)})
df <- data.frame(Reduce(cbind, normals))
colnames(df) <- 1:ncol(df)
df$mean <- rowMeans(df) # <- this is the average I am generating
df2 <- stack(df)
ggplot(df2, aes(values, fill = ind)) + geom_density(alpha = 0.000001) +theme(legend.position="none")

Por último, mi pregunta es cómo correctamente promediar los posteriors generados por mi modelo bajo diferentes semillas.

Answer 1

Basado en @ whuber comentario. La solución se encuentra promediando todo el vector de posteriors (de todas las ejecuciones) lo que me permite promediar el densidades no el variables aleatorias como hice en mi pregunta original.

library(ggplot2)
set.seed(10100101)
N <- 150;R <- 25
normals <- lapply(1:R,function(x){rnorm(n = N, mean = 0, sd   = 1)})
df <- data.frame(Reduce(cbind, normals))
colnames(df) <- 1:ncol(df)
df2 <- stack(df)

ggplot() +
  geom_density(data = df2, aes(x = values, fill = ind), alpha = 0) +
  geom_density(data = df2, aes(x = values), alpha = 0.3, colour="red",size=1.5)  + 
  theme(legend.position="none")