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 = [[1.0, 0.6, 0.9],
              [0.6, 1.0, 0.5],
              [0.9, 0.5, 1.0]]

cor_matrix = np.array(cor_matrix)
L = np.linalg.cholesky(cor_matrix)

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

[[ 1.          0.34450587  0.57515737]
 [ 0.34450587  1.          0.1488504 ]
 [ 0.57515737  0.1488504   1.        ]]