Am I wrong in assuming that the bayesian estimate has to converge to the true mean? I wrote a small simulation in R wherein I simulate two slightly different groups that are normally distributed. I then simulated a series of ten identical experiments where 25 new participants are sampled for each group for each experiment.
I start with a simple bayesian model with a noniformative prior for the effect of group on Y. Then:


*

*I take the posterior mean and sd of the model

*construct a new prior from the posterior mean and sd

*feed this prior into the model for the next experiment


and so on and so on. basically, the posterior of experiment t-1 becomes the prior for experiment t.
My intuition would suggest that over time, the estimate should slowly converge towards the true mean of the population. However, this is not what happens. Am I missing something?
Here is the R-Code in question, I added some comments:
#########################################################################################
# loading the libraries
#########################################################################################

if (!require(tidyverse)) install.packages('tidyverse')
#if (!require(cowplot)) install.packages('cowplot')
if (!require(brms)) install.packages('brms')
library(tidyverse)
#library(cowplot)
library(brms)



#########################################################################################
# Simulating Data
#########################################################################################

# difference between means of two groups
difference <- 0.2

# variation within group
variation <- 1

# number of experiments
experiments <- 10

# samples per experiment
n <- 25

#------------------------------------------------------------------------------------
#simulation of new data for the two groups
#------------------------------------------------------------------------------------

mysample <- c(rnorm(n,0,variation),rnorm(n,difference,variation))
df <- data.frame(mysample,group=rep(c(1,0),each=n))

#------------------------------------------------------------------------------------
#variables reserved for the for loop
#------------------------------------------------------------------------------------

mofs <- rep(NA,experiments)
mypmean <- rep(NA,experiments)
mypsd <- rep(NA,experiments)

# dimensions of the plot
minva <- difference-(variation*3)
maxva <- difference+(variation*3)


#########################################################################################
# Kickstart the Model
#########################################################################################



# We first need to kickstart the model with some rather broad prior:
bmodel1 <- brm(formula=mysample~group,
               prior=c(
                 set_prior("normal(0,50)",class="b",coef="group") #set an uninformative prior
               ),
               data = df,
               warmup = 200, iter = 1000, chains = 4,
               control = list(adapt_delta = 0.95),
               cores = getOption("mc.cores", 3L))

#########################################################################################
# Extract the posterior summary
#########################################################################################

# extract posterior
param <- posterior_summary(bmodel1)

#create a new Prior from the posterior
newprior <- paste0("normal(",param[[2]],",",param[[3]],")")


#########################################################################################
# Run the specified number of experiments
#########################################################################################

for (i in c(1:experiments)) {

# two groups with mean difference as specified and a variation as specified. combined as a data frame
  mysample <- c(rnorm(n,0,variation),rnorm(n,difference,variation))
  df <- data.frame(mysample,group=rep(c(1,0),each=n))

# the mean difference between the groups of the sample "mofs" is saved just to see how it varied in each experiment  
  mofs[i] <- mean(df$mysample[which(df$group==1)]-
                  df$mysample[which(df$group==0)])


#------------------------------------------------------------------------------------
#The kickstarted model us updated with new data AND a new prior (the posterior of the experiment before)
#------------------------------------------------------------------------------------
    tmp <- update(bmodel1,prior=
                  set_prior(newprior),newdata = df)

# save the estimated difference between both groups  and the sd
  mypmean[i] <- param[[2]]
  mypsd[i] <- param[[3]]

#------------------------------------------------------------------------------------
# create a new prior from the posterior
#------------------------------------------------------------------------------------
  param <- posterior_summary(tmp)
  newprior <- paste0("normal(",param[[2]],",",param[[3]],")")
}
#####################################################################################
# Plot the results
####################################################################################

df <- data.frame(estimatedDif=mypmean,
                 index=c(1:experiments),
                 sd=mypsd,
                 sampledDif=mofs)

x11()
ggplot(df,aes(index,estimatedDif))+
  geom_point(aes(index,sampledDif),size=3,alpha=0.3,color="lightblue")+
  geom_point(alpha=1)+
  ylim(c(minva,maxva))+
  geom_hline(aes(yintercept=difference*-1),color="green")+
  ggtitle(paste0("Sampled mean (blue) and Estimated (black) for ",experiments," Experiments and n=",n))

 A: Okay. So there where three things that I did wrong:


*

*I did not include a prior for the intercept and sigma

*I added the wrong parameter (sigma) to the prior for the group instead if its SE

*The "update" function needs to recompile the model to incoorporate the new priors, otherwise it will just proceed with the old uniformative prior. 


Of course, letting brms recompile the model for each subsequent update takes a lot more time so that 300 simulations are not feasible for my very, very old laptop. So I just proceeded with ten experiments, an effect of 0.2 n=25 and a variation of 1:

The changes in the code are
    tmp <- update(bmodel1,prior=c(
      set_prior(newpriorInt,class="Intercept"),
      set_prior(newpriorG,class="b",coef="group"),
      set_prior(newpriorSig,class="sigma")
    ),newdata = df)



#------------------------------------------------------------------------------------
# create a new prior from the posterior
#------------------------------------------------------------------------------------
  param <- posterior_summary(tmp)

  newpriorInt <- paste0("normal(",param[[1]],",",param[[1,2]],")")
  newpriorG <- paste0("normal(",param[[2,1]],",",param[[2,2]],")")
  newpriorSig <- paste0("normal(",param[[3,1]],",",param[[3,2]],")")

