Hi there, I would like to use brms to estimate mean and variance per group. I made it work in stan but would like to use brms to do the same thing. Here is the stan code:
data {
int<lower=1> Ntotal ;
int<lower=1> Ngroup ;
int<lower=1> Nsample ; // for priors
real y[Ntotal] ;
int group[Ntotal] ;
real meanY ;
real sdY ;
}
transformed data {
real normalSigma ;
normalSigma = sdY*10 ;
}
parameters {
real mu[Ngroup] ;
real<lower=0> sigma[Ngroup] ;
real mean_mu ;
real<lower=0> sigma_mu ;
}
model {
// hyper
mean_mu ~ normal(meanY, 100);
sigma_mu ~ uniform(1, 100) ;
// priors
sigma ~ cauchy( sdY , 100 ) ;
mu ~ normal( mean_mu, sigma_mu ) ;
for (i in 1:Ntotal){
y[i] ~ normal( mu[group[i]] , sigma[group[i]] ) ;
}
}
generated quantities{
vector[Ntotal] logLikelihood ; // for model comparison
vector[Ntotal] y_post ; // for posterior checking
vector[Nsample] mean_mu_prior ; // for checking priors
vector[Nsample] sigma_mu_prior ;
vector[Nsample] sigma_prior ;
vector[Nsample] mu_prior ;
for (i in 1:Ntotal){
logLikelihood[i] = normal_lpdf(y[i] | mu[group[i]], sigma[group[i]]) ;
y_post[i] = normal_rng(mu[group[i]], sigma[group[i]]) ;
}
for (i in 1:Nsample){
mean_mu_prior[i] = normal_rng(meanY, 100) ;
sigma_mu_prior[i] = uniform_rng(1, 100) ;
sigma_prior[i] = cauchy_rng(sdY, 100) ;
mu_prior[i] = normal_rng(mean_mu_prior[i], sigma_mu_prior[i]) ;
}
}
In brms: I have the code:
fit1 <- brm(bf(deltaZ ~ Overlap_ID -1, sigma ~ Overlap_ID-1 ) , data=residual_df, family = gaussian(), warmup=1000, iter=5000, chains =4, inits=“random”, cores=4)
Basically, I would like brms to estimate mu[i] and sigma[i] per group (Overlap_ID). But the estimates of the sigma’s are way too small than what I got from Stan.
Thanks.