# Generated MVN not what I expect

**URL:** <https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910>\
**Category:** General\
**Created:** [December 10, 2018, 11:58am UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910 "2018-12-10T11:58:16Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 11:58am UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/1 "2018-12-10T11:58:16Z")

</div>

I wrote a short piece of code to generate two variables from a multivariate normal but the covariances of the generated distributions are not what I expect. Hoping someone can help me to understand where I’m going wrong. Here is my stan code:

```stan
data{
    int<lower=0> N; // Total number of rows
}

generated quantities{
    matrix[2, 2] Sigma; // Square form of covariance matrix
    vector[2] zeroes[N];
    vector[2] y_ppc [N]; // model predictions for posterior predictive check
    
    // Specify covariance matrix
    Sigma[1,1] = 1;
    Sigma[2,1] = 0.05;
    Sigma[1,2] = 0.05;
    Sigma[2,2] = 0.15;

    // Generate y_ppc    
    for (n in 1:N) {
        zeroes[n] = rep_vector(0, 2);
    }

    y_ppc = multi_normal_rng(zeroes, Sigma);
}

```

Ok I run that through R and extract the generated variables:

```
N = 1000
datalist <- list(N)

fit1 <- stan(file="Stan models/gen_mvn.stan",
             data = datalist, chains = 3,
             seed=408, refresh=25, iter = 500 , algorithm = "Fixed_param"
) 
params = extract(fit1)
y <- data.frame(t(sapply(1:datalist$N, function(x) colMeans(params$y_ppc[, x, 1:2]) )))

```

Great. Runs without hitch. However when I look at the covariance of the generated variables:

```
> format( cov(y), digits = 2)
   X1 X2       
X1 "0.00145" "0.00007"
X2 "0.00007" "0.00019"

```

This is different to covariance matrix I specified in the stan code which was:

```
1.0 0.05
0.05 0.15

```

Umm… why is this different?

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 1:45pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/2 "2018-12-10T13:45:03Z")

</div>

Wrong calculation, try this instead (I don’t recall if the `byrow` is right or not, easy to check…):

```
library(rstan)
library(magrittr)
N = 1000
datalist <- list(N)

fit1 <- rstan::stan(file="mvr.stan",
             data = datalist, chains = 3,
             seed=408, refresh=25, iter = 500 , algorithm = "Fixed_param"
)
params = rstan::extract(fit1)

est_cov = matrix(
  data = apply(params$y_ppc, 1, cov) %>% apply(1, mean), 
  ncol=2, byrow=TRUE)

```

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 2:34pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/3 "2018-12-10T14:34:21Z")

</div>

Thanks @sakrejda. Now I’d best check I’ve been extracting my other parameters across the correct dimensions…eep

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 2:47pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/4 "2018-12-10T14:47:05Z")

</div>

yeah, don’t we all :/ also: transform before you summarize! :)

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 4:14pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/5 "2018-12-10T16:14:42Z")

</div>

Hey again. Sorry I’ve gotten myself confused now. I adapted your code to try extract the estimates:

```
y <- data.frame( matrix(
    data = apply(params$y_ppc1, 1, c) %>% apply(1, mean), 
    ncol=2, byrow=TRUE) )

```

but the distributions are again much tighter thasn i expect:

```
> summary(y)
       X1 X2            
 Min. :-0.116121 Min. :-1.040e-01  
 1st Qu.:-0.013724 1st Qu.:-1.340e-02  
 Median : 0.000992 Median :-1.292e-03  
 Mean : 0.001208 Mean :-3.142e-05  
 3rd Qu.: 0.016030 3rd Qu.: 1.289e-02  
 Max. : 0.121384 Max. : 1.104e-01

```

Perhaps I am more fundamentally confused in some way. For clarity what I actually want is the equivalent of this:

```
library(MASS)
M <- matrix(c(1, 0.05, 0.05, 0.15), nrow=2, ncol=2)
Data=mvrnorm(n = 800, mu = c(0, 0),Sigma = M,empirical = F)

> cov(Data)
           [,1] [,2]
[1,] 1.00669272 0.01918344
[2,] 0.01918344 0.14648667

> summary(Data)
       V1 V2           
 Min. :-4.80409 Min. :-1.153805  
 1st Qu.:-0.70604 1st Qu.:-0.285007  
 Median :-0.06849 Median :-0.008895  
 Mean :-0.05267 Mean :-0.014432  
 3rd Qu.: 0.65333 3rd Qu.: 0.227168  
 Max. : 2.81348 Max. : 1.202469  
> 

```

I’m sure I’m going wrong somehow but I don’t get where ?

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 7:11pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/6 "2018-12-10T19:11:04Z")

</div>

Your margins are: 1) iteration, 2) N, 3) 1:2

You want to apply the ‘cov’ function once for each iteration, that gets you one covariance matrix per iteration, then you want to summarise that (e.g.-using mean, or quantiles, or whatever

I don’t see why you write: `apply(params$y_ppc1, 1, c)`, that should be `cov` not `c`.

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 7:29pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/7 "2018-12-10T19:29:45Z")

</div>

> [@sakrejda](#):
>
> I don’t see why you write: `apply(params$y_ppc1, 1, c)` , that should be `cov` not `c` .

I’m trying to extract the two estimates y1, y2 so this is `c()` not `cov()`.  
My own code for this was:

```no-highlight
y <- data.frame(t(sapply(1:datalist$N, function(x) colMeans(params$y_ppc[, x, 1:2]) )))

```

> [@sakrejda](#):
>
> You want to apply the ‘cov’ function once for each iteration, that gets you one covariance matrix per iteration, then you want to summarise that (e.g.-using mean, or quantiles, or whatever

This is perhaps the crux of the issue. I actually don’t care about the covariance matrix per iteration. I care about the overall covariance for the mean extracted y vectors. This is is what the `mvrnorm` gives me. The Stan code seems to be doing something different - or i’m extracting things wrong ?

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 7:41pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/8 "2018-12-10T19:41:39Z")

</div>

With ‘c’ you’re losing the matrix structure so no covariance, maybe either cbind or rbind depending on which margin is which

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 8:33pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/9 "2018-12-10T20:33:52Z")

</div>

Arrrgghhh. I think I understand the problem. It’s not the way I’m summarising the data. The covariance structure holds per iteration. But averaged over all the iterations it reduces to the mean - i.e. zero, for any given y1 or y2. Ehhh, thats not really what I wanted.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 10:22pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/10 "2018-12-10T22:22:32Z")

</div>

… It should be smaller but not zero.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [December 10, 2018, 11:16pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/11 "2018-12-10T23:16:40Z")

</div>

Should’ve fired the second neuron before posting, you’re completely correct.

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [December 10, 2018, 11:24pm UTC](https://discourse.mc-stan.org/t/generated-mvn-not-what-i-expect/6910/12 "2018-12-10T23:24:26Z")

</div>

> [@sakrejda](#):
>
> Should’ve fired the second neuron before posting, you’re completely correct.

ha i know the feelin i’ve been puzzled by this for days 🤦‍♂️