在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()

结果:

另请参阅SciPy Cookbook中的CorrelatedRandomSamples.

Answer 2

如果Cholesky将协方差矩阵分解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分解

如果要生成两个系列,X并且Y具有特定(Pearson)相关系数(例如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.

例如,下面我们产生对相关的系列,X和Y1000次.然后我们绘制相关系数的直方图:

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.