scipy

Question

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

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

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

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

我尝试过的

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

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

看来Iman, Conover paper是标准方法

Answer 1

看来像Metropolis-Hastings算法这样的基于拒绝的抽样方法才是你想要的。 Scipy 可以用它的 scipy.optimize.basinhopping 函数实现这样的方法。

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

剩下的技巧就是：

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

鉴于 (1) 的结果，(2) 应该很简单。

Answer 2

如果您只想通过高斯 Copula (*) 进行关联，则可以使用 numpy 和 scipy.

分几步计算

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

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

Answer 3

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

Answer 4

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 作为函数列表传递，每个函数都有一个输入和一个 return 值。

在我的特定用例中，我需要生成多变量-t-distributed 随机变量（除了 normal-distributed 之外），请参阅此答案以了解如何执行此操作：。此外，我在 back-transform 部分使用了 scipy.stats.t.cdf。

在我的特定用例中，所需的分布是代表预期财务损失的经验分布。然后必须将最终数据点加在一起以获得所有的总财务损失 individual-but-correlated 金融事件。因此，np.array(...).T 实际上在我的代码库中被 sum(...) 取代。

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

scipy - generate random variables with correlations

python

numpy