在 Python (3.3) 中生成相关数据

Posted

技术标签:

【中文标题】在 Python (3.3) 中生成相关数据【英文标题】:Generate correlated data in Python (3.3) 【发布时间】:2013-04-08 03:04:21 【问题描述】:

在 R 中有一个函数(cm.rnorm.cor,来自包 CreditMetrics),它获取样本数量、变量数量和相关矩阵以创建相关数据。

Python 中有没有等价物?

【问题讨论】:

对不起,我的错,在 Python 3.3 上。 哇...最近真的添加了支持!谢谢你提醒我。 @Dualinity,以一种不那么幽默的语气,除了来自 Blender 的大量软件包之外,我真的建议您尝试 Python(X,Y)。它是一组用于科学开发的 Python 包 + IPython + Great IDE,称为 spyder。 code.google.com/p/pythonxy 哈,好吧,我正在使用 Emacs...所以我在“IDE”部门负责... 我不敢相信在尝试了这么久之后... Scipy 已经安装好了。这个问题仍然存在,但我会删除“没有 SciPy”的部分。 【参考方案1】:

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


如果您想生成两个系列,XY,带有特定的 (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。

例如,下面我们生成一对相关序列,XY,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。

【讨论】:

您知道如何让数据准确地具有 0.2 的相关性(容差很小)吗? numpy.random.multivariate_normal 在幕后做了什么?因为我将前者与 Cholesky 方法进行了比较,发现后者明显更快,特别是对于更大维度的数据(比如几千)。 cholesky 方法是否仅适用于某些特定类型的协方差矩阵?我的 cov 矩阵只是对角线,或者非常稀疏。 @Jason multivariate_normal 使用 SVD 而不是 Cholesky【参考方案2】:

numpy.randomGenerator类的方法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.
rng = np.random.default_rng()
y = rng.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.', alpha=0.25)
plt.plot(mu[0], mu[1], 'ro', ms=3.5)
plt.ylabel('y[1]')
plt.axis('equal')
plt.grid(True)

plt.subplot(2,2,3)
plt.plot(y[:,0], y[:,2], 'b.', alpha=0.25)
plt.plot(mu[0], mu[2], 'ro', ms=3.5)
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.', alpha=0.25)
plt.plot(mu[1], mu[2], 'ro', ms=3.5)
plt.xlabel('y[1]')
plt.axis('equal')
plt.grid(True)

plt.show()

结果:

另请参阅 SciPy Cookbook 中的 CorrelatedRandomSamples。

【讨论】:

以上是关于在 Python (3.3) 中生成相关数据的主要内容,如果未能解决你的问题,请参考以下文章

在python中生成具有三个类的3个圆形数据集

如何忽略 numpy 数组中的 NaN 数据点并在 Python 中生成规范化数据?

R在数据框中生成非重复对,避免相同的组成员

阻止Python在shebang中生成pyc文件[重复]

如何在 Python 中生成唯一 ID? [复制]

如何在 Python 中生成动态(参数化)单元测试?