Search code examples
rbayesianwinbugsstanrstan

Rstan code for simple multivariate linear model


I'm trying to use Rstan to fit an example model from Christensen, Johnson, Branscum, and Hanson's Bayesian Ideas and Data Analysis: An Introduction for Scientists and Statisticians. The authors use WinBUGS, so some adaptation is necessary. The data are here and the WinBUGS code is copied at the bottom of this post. This is a very simple model but I'm a complete beginner, and I can't figure out how to get around the error I'm getting. My Stan code is as follows:

data {
  int N_subjects;
  int N_items;
  matrix[N_subjects,N_items] y;
}

parameters {
  vector[N_items] mu;
  real<lower=0> sigma;
  real<lower=-1,upper=1> rho;
}

transformed parameters {
  cov_matrix[N_items] Sigma;
  for (j in 1:N_items)
    for (k in 1:N_items)
      Sigma[j,k] <- pow(sigma,2)*pow(rho,step(abs(j-k)-0.5));
}

model {
  sigma ~ uniform(0,100);
  rho ~ uniform(0,1);
  mu ~ multi_normal(0,100);
  for (i in 1:N_subjects)
    y[i] ~ multi_normal(mu,Sigma);
}

The parser throws the following error:

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'model' with error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
no matches for function name="multi_normal_log"
    arg 0 type=vector
    arg 1 type=int
    arg 2 type=int
available function signatures for multi_normal_log:
0.  multi_normal_log(vector, vector, matrix) : real
1.  multi_normal_log(vector, row vector, matrix) : real
2.  multi_normal_log(row vector, vector, matrix) : real
3.  multi_normal_log(row vector, row vector, matrix) : real
4.  multi_normal_log(vector, vector[1], matrix) : real
5.  multi_normal_log(vector, row vector[1], matrix) : real
6.  multi_normal_log(row vector, vector[1], matrix) : real
7.  multi_normal_log(row vector, row vector[1], matrix) : real
8.  multi_normal_log(vector[1], vector, matrix) : real
9.  multi_normal_log(vector[1], row vector, matrix) : real
10.  multi_normal_log(row vector[1], vector, matrix) : real
11.  multi_normal_log(row vector[1], row vector, matrix) : real
12.  multi_normal_log(vector[1], vector

(I think) I understand that the parser is telling me I'm trying to pass inappropriate data types to the multi_normal function in the model block, but I can't figure out where this is originating from. I suspect I'm doing something wrong in defining the covariance matrix, but it appears that more than one argument has an incorrect data type...

WinBUGS code I'm modeling my Stan code on:

model{
for(i in 1:30){
for(j in 1:6){
logy[i,j] <- log(y[i,j])
}
}
for(i in 1:30){logy[i,1:6]~dmnorm(m[1:6],precision[1:6,1:6])}
for(j in 1:6){
for(k in 1:6){
covariance[j,k] <- sigma2*pow(rho, step(abs(j-k)-0.5))
}
}
for(i in 1:6){ m[i] <- mu }
precision[1:6,1:6] <- inverse(covariance[1:6,1:6])
sigma ~ dunif(0,100)
mu ~ dnorm(0,0.001)
L <- -1/(6-1)
rho ~ dunif(L,1)
sigma2 <- sigma*sigma
tau <- 1/sigma2
}

Solution

  • The error comes from

    mu ~ multi_normal(0,100);
    

    as you are passing a vector mu, integer 0, and integer 100. I suppose you want

    mu ~ normal(0,100);
    

    which treats the elements of mu as independent and identically normally distributed with mean 0 and standard deviation 100.