How to use the Cholesky decomposition, or an alternative, for correlated data simulation I use Cholesky decomposition to simulate correlated random variables given a correlation matrix. The thing is, the result never reproduces the correlation structure as it is given. Here is a small example in Python to illustrate the situation.
import numpy as np    

n_obs = 10000
means = [1, 2, 3]
sds = [1, 2, 3] # standard deviations 

# generating random independent variables 
observations = np.vstack([np.random.normal(loc=mean, scale=sd, size=n_obs)
                   for mean, sd in zip(means, sds)])  # observations, a row per variable

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

print(np.corrcoef(L.dot(observations))) 

This prints:
[[ 1.          0.34450587  0.57515737]
 [ 0.34450587  1.          0.1488504 ]
 [ 0.57515737  0.1488504   1.        ]]

As you can see, the post-hoc estimated correlation matrix drastically differs from the prior one. Is there a bug in my code, or is there some alternative to using the Cholesky decomposition? 
Edit
I beg your pardon for this mess. I didn't think there was an error in the code and/or in the way Cholesky decomposition was applied due to some misunderstanding of the material I had studied before. In fact I was sure that the method itself was not meant to be precise and I had been okay with that up until the situation that made me post this question. Thank you for pointing at the misconception I had. I've edited the title to better reflect the real situation as proposed by @Silverfish.
 A: People would be likely find your error much faster if you explained what you did with words and algebra rather than code (or at least writing it using pseudocode). 
You appear to be doing the equivalent of this (though possibly transposed): 


*

*Generate an $n\times k$ matrix of standard normals, $Z$

*multiply the columns by $\sigma_i$ and add $\mu_i$ to get nonstandard normals

*calculate $Y=LX$ to get correlated normals.
where $L$ is the left Cholesky factor of your correlation matrix.
What you should do is this:


*

*Generate an $n\times k$ matrix of standard normals, $Z$

*calculate $X=LZ$ to get correlated normals.

*multiply the columns by $\sigma_i$ and add $\mu_i$ to get nonstandard normals
There's many explanations of this algorithm on site. e.g. 
How to generate correlated random numbers (given means, variances and degree of correlation)?
Can I use the Cholesky-method for generating correlated random variables with given mean?
This one discusses it directly in terms of the desired covariance matrix, and also gives an algorithm for getting a desired sample covariance:
Generating data with a given sample covariance matrix
A: The approach based on the Cholesky decomposition should work, it is described here
and is shown in the answer by 
Mark L. Stone posted almost at the same time that this answer.
Nevertheless, I have sometimes generated draws from the multivariate Normal distribution 
$N(\vec\mu, \Sigma)$ as follows:
$$
Y = Q X + \vec\mu \,, \quad \hbox{with}\quad Q=\Lambda^{1/2}\Phi \,,
$$
where $Y$ are the final draws, $X$ are draws from the univariate standard Normal distribution, $\Phi$ is a matrix containing the normalized eigenvectors of the target matrix $\Sigma$ and $\Lambda$ is a diagonal matrix containing the eigenvalues of $\Sigma$ arranged in the same order as the eigenvectors in the columns of $\Phi$.
Example in R (sorry I'm not using the same software you used in the question):
n <- 10000
corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
set.seed(123)
SigmaEV <- eigen(corM)
eps <- rnorm(n * ncol(SigmaEV$vectors))
Meps <- matrix(eps, ncol = n, byrow = TRUE)    
Meps <- SigmaEV$vectors %*% diag(sqrt(SigmaEV$values)) %*% Meps
Meps <- t(Meps)
# target correlation matrix
corM
#      [,1] [,2] [,3]
# [1,]  1.0  0.6  0.9
# [2,]  0.6  1.0  0.5
# [3,]  0.9  0.5  1.0
# correlation matrix for simulated data
cor(Meps)
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.6002078 0.8994329
# [2,] 0.6002078 1.0000000 0.5006346
# [3,] 0.8994329 0.5006346 1.0000000


You may be also interested in 
this post
and this post.
A: As others have already shown: cholesky works. Here a piece of code which is very short and very near to pseudocode: a codepiece in MatMate:
Co = {{1.0, 0.6, 0.9},  _
      {0.6, 1.0, 0.5},  _
      {0.9, 0.5, 1.0}}           // make correlation matrix


chol = cholesky(co)              // do cholesky-decomposition           
data = chol * unkorrzl(randomn(3,100,0,1))  
                                 // dot-multiply cholesky with random-
                                 // vectors with mean=0, sdev=1  
                                 //(refined by a "decorrelation" 
                                 //to remove spurious/random correlations)   


chk = data *' /100               // check the correlation of the data
list chk

1.0000  0.6000  0.9000
0.6000  1.0000  0.5000
0.9000  0.5000  1.0000

A: Python code:
import numpy as np

# desired correlation matrix
cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

# build some signals that will result in the desired correlation matrix
X = L.dot(np.random.normal(0,1,(3,1000))) # the more the sample (1000) the better

# estimate their correlation matrix
np.corrcoef(X)
array([[1.        , 0.58773667, 0.8978625 ],
       [0.58773667, 1.        , 0.47424997],
       [0.8978625 , 0.47424997, 1.        ]])

# Very good approxiamation :)


