Scipy - 生成具有相关性的随机变量

Question

Scipy - 生成具有相关性的随机变量

pythonnumpyscipy

16

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

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

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

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

我的尝试

这个答案和这个答案似乎表明“copulas”是一个答案，但我在scipy中没有看到任何参考。

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

似乎 Iman，Conover的论文是标准方法。

- MikeRand

这是您要找的吗？https://dev59.com/m2Qo5IYBdhLWcg3wfPQa#16025584 - unutbu

1

适用于普通变量...我有其他分布。 - MikeRand

看起来推荐的方法（Iman-Conover）使用多元正态分布来实现我所需要的功能，因此我认为您的评论可能是最终解决方案的重要部分（这可能是我需要手动构建的东西）。 - MikeRand

你能否分享一下你编写的用于生成具有相关性的随机变量的Python代码？ - Ark

你的问题不完整，因为边际未指定。根据Sklar定理，分布函数由其边际分布和其Copula完全确定。各种Copula将产生相关性：尽管高斯Copula是一种特定的选择，但还有许多其他选择。 - Michael Baudin

n个变量的集合不就是边缘分布的集合吗？ - MikeRand

4个回答

2

似乎像马尔科夫链蒙特卡罗采样方法（如Metropolis-Hastings算法）是您想要的。Scipy可以通过其scipy.optimize.basinhopping函数实现这些方法。

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

然后，剩下的技巧将是：

将文本翻译成中文：

找出具有所需边缘分布形式的联合N维概率密度函数的形式，但是具有所需的相关矩阵。这对于高斯分布来说很容易做到，其中所需的相关矩阵和均值向量就足以定义分布。如果您的边际具有简单的表达式，您可能可以通过一些直接但繁琐的代数找到这个pdf。此论文引用了几篇其他的论文，这些论文都在讨论您所说的问题，我相信还有更多。
制定一个函数用于最小化basinhopping，使得它被接受的“最小值”等于您定义的这个pdf的样本。

根据（1）的结果，（2）应该很容易做到。

- stochastic

1

如果您已经有一个正半定相关矩阵 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    ]

- Jean A.

1

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

0

import typing

import numpy as np
import scipy.stats


def run_gaussian_copula_simulation_and_get_samples(
        ppfs: typing.List[typing.Callable[[np.ndarray], np.ndarray]],  # List of $num_dims percentile point functions
        cov_matrix: np.ndarray,  # covariance matrix, shape($num_dims, $num_dims)
        num_samples: int,  # number of random samples to draw
) -> np.ndarray:
    num_dims = len(ppfs)

    # Draw random samples from multidimensional normal distribution -> shape($num_samples, $num_dims)
    ran = np.random.multivariate_normal(np.zeros(num_dims), cov_matrix, (num_samples,), check_valid="raise")

    # Transform back into a uniform distribution, i.e. the space [0,1]^$num_dims
    U = scipy.stats.norm.cdf(ran)

    # Apply ppf to transform samples into the desired distribution
    # Each row of the returned array will represent one random sample -> access with a[i]
    return np.array([ppfs[i](U[:, i]) for i in range(num_dims)]).T  # shape($num_samples, $num_dims)

# Example 1. Uncorrelated data, i.e. both distributions are independent
f1 = run_gaussian_copula_simulation_and_get_samples(
    [lambda x: scipy.stats.norm.ppf(x, loc=100, scale=15), scipy.stats.norm.ppf],
    [[1, 0], [0, 1]],
    6
)
# Example 2. Completely correlated data, i.e. both percentiles match
f2 = run_gaussian_copula_simulation_and_get_samples(
    [lambda x: scipy.stats.norm.ppf(x, loc=100, scale=15), scipy.stats.norm.ppf],
    [[1, 1], [1, 1]],
    6
)
np.set_printoptions(suppress=True)  # suppress scientific notation
print(f1)
print(f2)

关于这个函数的一些说明。np.random.multivariate_normal为我们做了很多繁重的工作，特别要注意的是我们不需要分解相关矩阵。ppfs被传递为一个函数列表，每个函数都有一个输入和一个返回值。

在我的特定用例中，我需要生成多元t-分布随机变量（除了正态分布），请参考此答案以了解如何实现：https://dev59.com/957ha4cB1Zd3GeqPnbQo#41967819。此外，我使用scipy.stats.t.cdf进行反变换部分。

在我的特定用例中，所需的分布是代表预期财务损失的经验分布。然后，最终数据点必须相加以获得所有单个但相关的财务事件的总财务损失。因此，在我的代码库中，np.array(...).T实际上被sum(...)替换。

- xjcl

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Josef · Accepted Answer

如果您只是想通过高斯Copula（*）进行相关性计算，则可以使用numpy和scipy在几个步骤中完成。

使用所需协方差创建多元随机变量，使用numpy.random.multivariate_normal函数并创建一个(nobs, k_variables)的数组。
对于每个列/变量，应用scipy.stats.norm.cdf将正态转换为均匀分布，以获得均匀边际分布。
应用dist.ppf将均匀分布转换为所需的分布，其中dist可以是scipy.stats中的分布之一。

（*）高斯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做过这个，但现在没有任何脚本或函数。）