# Multi_normal_rng covariance matrix unexpected behaviour

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

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.