for the case of normally-distributed data, the process is very simple:

(1) decide on the correlation matrix that you want

(2) do some matrix decomposition on it (cholesky/PCA i think are the most popular ones)

(3) multiply the decomposition of the correlation matrix times the matrix of (normally-distributed) variables

the result is a nice multivariate normal distribution with the correlation matrix you intended in (1).

i have found a couple of articles now that apply the same algorithm but for non-normal data. i reproduced their examples in R and i can see they work as far as preserving the correlation structure goes but that there are important changes in the moments of the distribution. now, i already knew that the matrix multiplication would alter the moments of the data... but what i didn't suspect is that the new datasets look...well... for lack of a better word, more "normally distributed".

lemme show you what i mean:

Code:

```
set.seed(123)
R <- matrix(c(1.0, 0.5, 0.5,
0.5, 1.0, 0.5,
0.5, 0.5, 1.0),3,3)
n<- 5000
unf <- matrix(c(runif(n), runif(n), runif(n)),n,3)
datz <- unf%*%chol(R)
```

**BEFORE THE TRANSFORMATION**

**AFTER THE TRANSFORMATION**

and i can even just look at the descriptives and see just how much it changed:

Code:

```
library(psych)
> describe(unf)
var n mean sd median trimmed mad min max range skew kurtosis se
1 1 5000 0.5 0.29 0.50 0.5 0.37 0 1 1 0.01 -1.20 0
2 2 5000 0.5 0.29 0.51 0.5 0.38 0 1 1 -0.03 -1.22 0
3 3 5000 0.5 0.29 0.50 0.5 0.37 0 1 1 0.00 -1.21 0
> describe(datz)
var n mean sd median trimmed mad min max range skew kurtosis se
1 1 5000 0.50 0.29 0.50 0.50 0.37 0.00 1.00 1.00 0.01 -1.20 0
2 2 5000 0.68 0.29 0.69 0.68 0.33 0.00 1.35 1.35 -0.02 -0.80 0
3 3 5000 0.80 0.29 0.81 0.80 0.32 0.06 1.53 1.48 0.00 -0.65 0
```

i mean, variable #3 reduced its kurtosis from -1.21 to -0.65! that's almost less than half the kurtosis!

i also explored the beta distribution with parameters shape1=shape2=0.5 as an example of another symmetric distirbution and the same thing appeared:

**BEFORE THE TRANSFORMATION:**

**AFTER THE TRANSFORMATION:**

i mean, here it was so extreme it even lost its bimodality!

i tried it with skewed distributions as well. i can see they don't lose their skewness, but the skewness is also considerably reduced after they get transformed.

i've been looking around all day on the internet but it doesn't seem like anyone has explored this method within the context of non-normal data. more importantly, nobody seems to explain why what initially looks like non-normal distributions become more "bell-shaped" after the transformation, ESPECIALLY if they are symmetric distributions.

ideas?