在Python中生成相关数据（3.3）

Question

在R中有一个函数（cm.rnorm.cor，来自程序包CreditMetrics），该函数用于创建相关数据的样本数量，变量数量和相关矩阵。

Python中是否有等价物？

Answer 1

numpy.random.multivariate_normal是你想要的功能。

实施例：

import numpy as np 
import matplotlib.pyplot as plt 


num_samples = 400 

# The desired mean values of the sample. 
mu = np.array([5.0, 0.0, 10.0]) 

# The desired covariance matrix. 
r = np.array([ 
     [ 3.40, -2.75, -2.00], 
     [ -2.75, 5.50, 1.50], 
     [ -2.00, 1.50, 1.25] 
    ]) 

# Generate the random samples. 
y = np.random.multivariate_normal(mu, r, size=num_samples) 


# Plot various projections of the samples. 
plt.subplot(2,2,1) 
plt.plot(y[:,0], y[:,1], 'b.') 
plt.plot(mu[0], mu[1], 'ro') 
plt.ylabel('y[1]') 
plt.axis('equal') 
plt.grid(True) 

plt.subplot(2,2,3) 
plt.plot(y[:,0], y[:,2], 'b.') 
plt.plot(mu[0], mu[2], 'ro') 
plt.xlabel('y[0]') 
plt.ylabel('y[2]') 
plt.axis('equal') 
plt.grid(True) 

plt.subplot(2,2,4) 
plt.plot(y[:,1], y[:,2], 'b.') 
plt.plot(mu[1], mu[2], 'ro') 
plt.xlabel('y[1]') 
plt.axis('equal') 
plt.grid(True) 

plt.show() 

结果：

参见CorrelatedRandomSamples在SciPy的食谱。

Answer 2

如果乔莱斯基-分解协方差矩阵C为L L^T，并产生 独立随机向量x，然后Lx将与协方差 C随机向量。

import numpy as np 
import matplotlib.pyplot as plt 
linalg = np.linalg 
np.random.seed(1) 

num_samples = 1000 
num_variables = 2 
cov = [[0.3, 0.2], [0.2, 0.2]] 

L = linalg.cholesky(cov) 
# print(L.shape) 
# (2, 2) 
uncorrelated = np.random.standard_normal((num_variables, num_samples)) 
mean = [1, 1] 
correlated = np.dot(L, uncorrelated) + np.array(mean).reshape(2, 1) 
# print(correlated.shape) 
# (2, 1000) 
plt.scatter(correlated[0, :], correlated[1, :], c='green') 
plt.show() 

参考：见Cholesky decomposition

---

如果你想生成两个系列，X和Y，特别(Pearson) correlation coefficient（如0.2）：

rho = cov(X,Y)/sqrt(var(X)*var(Y)) 

你可以选择的协方差矩阵是

cov = [[1, 0.2], 
     [0.2, 1]] 

这使得cov(X,Y) = 0.2，和方差，var(X)和var(Y)都等于1，所以rho将等于0.2。

例如，下面我们生成1000对相关系列，X和Y,1000次。然后我们绘制的相关系数的柱状图：

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.stats as stats 
linalg = np.linalg 
np.random.seed(1) 

num_samples = 1000 
num_variables = 2 
cov = [[1.0, 0.2], [0.2, 1.0]] 

L = linalg.cholesky(cov) 

rhos = [] 
for i in range(1000): 
    uncorrelated = np.random.standard_normal((num_variables, num_samples)) 
    correlated = np.dot(L, uncorrelated) 
    X, Y = correlated 
    rho, pval = stats.pearsonr(X, Y) 
    rhos.append(rho) 

plt.hist(rhos) 
plt.show() 

正如你所看到的，相关系数一般都接近0.2，但对于任何给定的样本，该相关性将最有可能不是0.2完全相同。