## Generating a multivariate gaussian distribution using RcppArmadillo

Ahmadou Dicko — written Mar 12, 2013 — source

There are many ways to simulate a multivariate gaussian distribution assuming that you can simulate from independent univariate normal distributions. One of the most popular method is based on the Cholesky decomposition. Let’s see how `Rcpp` and `Armadillo` perform on this task.

The easiest way to perform a Cholesky distribution in R is to use the `chol` function in `R` which interface some fast `LAPACK` routines.

We will also use `MASS:mvrnorm` which implements it differently:

```Loading required package: MASS
```
```        [,1]    [,2]    [,3]
[1,]  1.0000  0.8851 -0.3830
[2,]  0.8851  1.0000 -0.4675
[3,] -0.3830 -0.4675  1.0000
```
```        [,1]    [,2]    [,3]
[1,]  1.0000  0.9106 -0.3016
[2,]  0.9106  1.0000 -0.4599
[3,] -0.3016 -0.4599  1.0000
```
```       [,1]    [,2]    [,3]
[1,]  1.000  0.9020 -0.3530
[2,]  0.902  1.0000 -0.4889
[3,] -0.353 -0.4889  1.0000
```
```[1]  9.850  4.911 -2.902
```
```[1] 10.051  5.046 -2.914
```
```[1]  9.825  4.854 -2.873
```

Now, let’s benchmark the different versions:

```Loading required package: rbenchmark
```
```                             test replications relative elapsed
3   mvrnormArma(10000, mu, sigma)          100    1.000   0.219
1      mvrnormR(10000, mu, sigma)          100    1.913   0.419
2 MASS::mvrnorm(10000, mu, sigma)          100    2.046   0.448
```

The `RcppArmadillo` function outperforms the MASS implementation and the naive R code, but more surprisinugly `mvrnormR` is slightly faster than `mvrnorm` in this benchmark.

To be fair, while digging into the `MASS::mvrnorm` code it appears that there are few code sanity checks ( such as the positive definiteness of `Sigma` ).