scipy - 生成具有相关性的随机变量答案

【问题标题】：scipy - generate random variables with correlationsscipy - 生成具有相关性的随机变量
【发布时间】：2015-02-27 22:26:34
【问题描述】：

我正在用 Python 实现一个基本的 Monte Carlo 模拟器，用于我正在尝试做的一些项目管理风险建模（基本上是 Crystal Ball / @Risk，但在 Python 中）。

我有一组n 随机变量（所有scipy.stats 实例）。我知道我可以使用rv.rvs(size=k) 从每个n 变量中生成k 独立观察结果。

我想通过指定一个n x n 半正定相关矩阵来引入变量之间的相关性。

在 scipy 中有没有干净的方法来做到这一点？

我的尝试

This answer 和 this answer 似乎表明“copulas”将是一个答案，但我在 scipy 中没有看到对它们的任何引用。

This link 似乎实现了我正在寻找的东西，但我不确定 scipy 是否已经实现了这个功能。我也希望它适用于非正态变量。

看来Iman, Conover paper是标准方法。

【问题讨论】：

这是您要找的吗？ stackoverflow.com/a/16025584/190597
适用于正态变量...我还有其他分布。
看来推荐的方法（Iman-Conover）使用多元法线来做我正在寻找的东西，所以我认为您的评论可能是最终解决方案的很大一部分（这可能是我必须手工构建的东西）。
您是否有机会分享您为生成具有相关性的随机变量而开发的 Python 代码？
您的问题不完整，因为未指定边缘。根据 Sklar 定理，分布函数完全由其边际分布及其 copula 指定。各种 copula 会产生相关性：虽然 Gaussian copula 是一个特定的选择，但还有很多其他选择。

标签： python numpy scipy

【解决方案1】：

您想要的是一种基于拒绝的抽样方法，例如 Metropolis-Hastings 算法。 Scipy 可以通过其scipy.optimize.basinhopping 函数实现此类方法。

基于拒绝的抽样方法允许您从任何给定的概率分布中抽取样本。这个想法是，您从另一个易于采样的“提案”pdf 中抽取随机样本（例如均匀分布或高斯分布），然后使用随机测试来确定提案分布中的该样本是否应该被“接受”为代表所需分布的样本。

剩下的技巧将是：

找出联合 N 维概率密度函数的形式，该函数在每个维度上具有所需形式的边际，但具有所需的相关矩阵。这对于高斯分布很容易做到，其中所需的相关矩阵和均值向量就是定义分布所需的全部。如果你的边际有一个简单的表达式，你可能会发现这个 pdf 有一些简单但乏味的代数。 This 论文引用了其他几个你正在谈论的内容，我敢肯定还有更多。
为basinhopping 制定一个函数，以使其最小化，使其接受的“最小”数量为您定义的此pdf 的样本。

鉴于 (1) 的结果，(2) 应该是直截了当的。

【讨论】：

【解决方案2】：

如果您只想通过高斯 Copula (*) 进行相关性，则可以使用 numpy 和 scipy 分几步计算。

创建具有所需协方差的多元随机变量numpy.random.multivariate_normal，并创建一个（nobs by k_variables）数组
应用scipy.stats.norm.cdf 将正态随机变量转换为均匀随机变量，为每列/变量获得均匀的边际分布
应用dist.ppf 将统一边距转换为所需的分布，其中dist 可以是scipy.stats 中的分布之一

(*) Gaussian copula 只是一种选择，当我们对尾部行为感兴趣时它不是最好的，但它是最容易使用的例如http://archive.wired.com/techbiz/it/magazine/17-03/wp_quant?currentPage=all

两个参考

https://stats.stackexchange.com/questions/37424/how-to-simulate-from-a-gaussian-copula

http://www.mathworks.com/products/demos/statistics/copulademo.html

（我之前可能在 python 中做过这个，但现在没有任何脚本或函数。）

【讨论】：

您知道任何类似的、内存效率更高的解决方案吗？我正在使用 'cov_matrix = toeplitz(rho**arange(p))' 执行此操作，但是当我达到高维度时会遇到内存错误。
如何在 python 中获得均匀的边缘分布？
@Ark 要获得均匀的边际分布，您可以跳过最后一步。
“创建（k_variables 的nobs）数组”是什么意思？ @约瑟夫
nobs 是观察次数。您创建每个维度 k_variables 的 nobs 随机变量。对于数据分析，我们通常在行中有观察结果，在列中有数据系列或变量。

【解决方案3】：

如果您已经有一个半正定相关矩阵 R [n x n]，则可以轻松构建一个以 R 作为输入的 NormalCopula。我将向您展示一个 n = 3 的示例。代码基于OpenTURNS library。

import openturns as ot

# you can replace this part by your matrix
dim = 3
R = ot.CorrelationMatrix (dim)
R[0,1] = 0.25
R[0,2] = 0.6
R[1,2] = 0.9

copula = ot.NormalCopula(R)

如果你想得到一个大小的样本，就写

size = 5
print(copula.getSample(size))
>>>    [ X0       X1       X2       ]
0 : [ 0.355353 0.76205  0.632379 ]
1 : [ 0.902567 0.984443 0.989552 ]
2 : [ 0.423219 0.811016 0.754304 ]
3 : [ 0.303776 0.471557 0.450188 ]
4 : [ 0.746168 0.918729 0.891347 ]

编辑 - 遵循@Michael_Baudin 的评论

当然，如果您想将边际分布设置为例如Beta 和 LogNormal 边际，它也是可能的：

X0 = ot.LogNormal(0.1, 1, 0)
X1 = ot.Beta()
X2 = ot.Uniform(1.0, 2.0)
distribution = ot.ComposedDistribution([X0,X1,X2], Original_copula)
print(distribution.getSample(size))
>>> [ X0         X1         X2         ]
0 : [  3.97678    0.158823   1.75635   ]
1 : [  1.18929   -0.554092   1.18952   ]
2 : [  2.59542    0.0751359  1.68599   ]
3 : [  1.33363   -0.18407    1.42241   ]
4 : [  1.34084    0.198019   1.6553    ]

【讨论】：

我建议扩展脚本并设置边缘分布，例如Beta 和 LogNormal 边际，因为问题提到“我也希望它适用于非正态变量。”。