How to build a simple hierarchical model with different samples in Stan/rstan?

Modeling
,
1

I’m a newbie in stan.

I’m trying to build a simple hierarchical model.

First of all I create five normally distributed samples for modelling with random mean and standard deviation.

library(rstan)

set.seed(10)

mus    = rnorm(5,2,2)
sigmas = rlnorm(5,3,1)

N=40
J=5
y = matrix(NA,nrow = N,ncol = J)

# make data matrix
for (i in 1:5) {
  y[,i] = rnorm(N,mus[i],sigmas[i])
}

model = stan_model(paste0(dir,"NormalHierarchical.stan"))
fit = sampling(model,list(N=N,J=J,y=y),iter=200,chains=1)

In the data level each sample has a normal distribution N(\mu,\sigma).
In the process level, \mu \sim N(\tau,\omega) and \sigma \sim logN(\lambda,\delta).

Where logN is the lognormal distribution.

And the top level priors are:
\tau \sim N(0,1)
\omega \sim logN(0,1)
\lambda \sim N(0,1)
\delta \sim logN(0,1)

My stan code is as follows and doesn’t work

data {
  int<lower=1> N; //salple length
  int<lower=1> J; //number of samples
  matrix[N,J] y;  //data matrix
}

parameters {
  vector[J] mu;
  vector[J] sigma;
  
  real tau; 
  real<lower=0> omega;
  real lambda; 
  real<lower=0> delta;
}

model {
  //priors
  tau ~ normal(0,1);
  omega ~ lognormal(0,1);
  lambda ~ normal(0,1);
  delta ~ lognormal(0,1);
  
  for(j in 1:J) {
    mu[j] ~ normal(tau,omega);
    sigma[j] ~ lognormal(lambda,delta);
    y[,j] ~ normal(mu,sigma);
  }
}

The error code

SAMPLING FOR MODEL 'NormalHierarchical' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: normal_lpdf: Location parameter has dimension = 5, expecting dimension = 40; a function was called with arguments of different scalar, array, vector, or matrix types, and they were not consistently sized;  all arguments must be scalars or multidimensional values of the same shape.  (in 'model31c0139f25ab_NormalHierarchical' at line 27)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                                                                                                                                                                                                              
[2] "  Exception: normal_lpdf: Location parameter has dimension = 5, expecting dimension = 40; a function was called with arguments of different scalar, array, vector, or matrix types, and they were not consistently sized;  all arguments must be scalars or multidimensional values of the same shape.  (in 'model31c0139f25ab_NormalHierarchical' at line 27)"
[3] "In addition: Warning message:"                                                                                                                                                                                                                                                                                                                                 
[4] "In readLines(file, warn = TRUE) :"                                                                                                                                                                                                                                                                                                                             
[5] "  incomplete final line found on '/.../NormalHierarchical.stan'"                                                                                                                                                                                                                                               
error occurred during calling the sampler; sampling not done

I’m trying to make the output for a vector of posteriors for each sample i.e. mu[1],mu[2],mu[3]… and sigma[1],sigma[2],sigma[3]…

Maybe a data matrix is not the right approach? As said, I’m only using stan since a few weeks.

2

The error indicates that

    y[,j] ~ normal(mu,sigma);

is expecting y[,] to be a vector with 5 elements. Changing

  matrix[N,J] y;  //data matrix

to

  matrix[J,N] y;  //data matrix

should fix it.