How to average several posteriors distributions from a Monte Carlo Simulation Say you produce several posteriors distributions from different runs of the same model under different seeds. That is to say you have something like the following:
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")


I would like to add to the plot the "average posterior" produced by the different runs of my model. At first I thought that I would be as easy as to apply the mean of the posteriors. However, I am getting something that converges to zero because I am simply averaging values that came from a standard normal distribution as you can see in the picture below:
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")


Finally, my questions is how to correctly average the posteriors generated by my model under different seeds.

Based on @whuber comment. The solution is found by averaging the entire vector of posteriors (of all runs) which allows me to average the densities not the random variables as I did in my original question.
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") 


