Multi_normal_rng covariance matrix unexpected behaviour

General
1

I get a strange behaviour from multi_normal_rng function.
In R, I do this:

library(MASS)

# simulate from a multinormal distribution
mu <- rep(0, 6)
S  <- matrix(c(2.99953664, 1.166298844, 1.001939952, 0, 1.164095096, 
  -0.69732028, 1.166298844, 3.12289268, -1.0965728, -0.999813992, 
  1.146612724, 0, 1.001939952, -1.0965728, 3.5697292, 1.154118408, 
  0, -0.89103026, 0, -0.999813992, 1.154118408, 3.42587208, 1.336106332, 
  1.229933084, 1.164095096, 1.146612724, 0, 1.336106332, 3.51555988, 
  1.21955636, -0.69732028, 0, -0.89103026, 1.229933084, 1.21955636, 
  2.78408564), 6, 6)

mvrnorm(100, mu, S) %>% apply(., 2, sd)
# I obtain the following standard deviation, so far so good:
# 1.648693 1.735414 2.027764 1.993323 1.999680 1.495854

Now I create the Stan code again to simulate from a multinormal distribution using the same parameters:


data {
  int<lower=0> N;
  int<lower=0> M;
  vector[M] mu;
  matrix[M,M] S;
}

generated quantities {
  array[N]vector[M] r_hat;
  for(i in 1:N) {
      r_hat[i] = multi_normal_rng(mu, S);
  }  
}

I compile, run and extract the r_hat values:.

fit <- mod$sample(data = list(N=100, M=6, mu=mu, S=S), chains = 1, iter_warmup = 250, iter_sampling = 250,fixed_param = TRUE)
fit$draws(variables = "r_hat") %>% colMeans() %>% matrix(., 100, 6) %>% apply(., 2, sd)
# However, I get the following std:
# [1] 0.1140576 0.1059731 0.1190921 0.1180816 0.1150848 0.1085626

I expected mvrnorm and multi_normal_rng to simulate the same values, with the same standard deviations.
What am I doing wrong? The way I extract the r_hat values?

thank you

2

You took the column means before taking the standard deviations. Since sd() is a nonlinear operation, you need to compute sd() iteration-wise and then summarize that distribution however you see fit, e.g. by taking its mean.