使用 SciPy 的 QMC (Sobel) 生成多个非同分布的准随机变量

Question

我在使用 SciPy 的准蒙特卡罗特征进行函数估计时遇到问题。我想为多个变量的函数构建生存函数，所以我使用这样的东西：

Y = c₁f₁(X₁) + c₂f₂(X₂) + c₃f₃(X₃)…

地点：

• c₁、c₂、c₃ 是常数

• f₁、f₂、f₃…是函数（大部分是线性的，但并不总是）

• X₁、X₂、X₃...是独立的，但不是同分布的随机变量。事实上，分布大多不同。

对于“正常”蒙特卡罗，我已经解决了问题，我只是生成具有不同分布的随机数，然后将它们放入上面的方程中。但对于准蒙特卡罗，我陷入了困境。

如果随机变量都是相同的分布，我可以这样做（x_rvs是我的随机变量，并且它们没有_dims）：

engine = scipy.stats.qmc.Sobol(d=no_dims, scramble=True)
dist = scipy.stats.fisk(3.9)
base = scipy.stats.sampling.NumericalInverseHermite(dist)
x_rvs = base.qrvs(size=size, qmc_engine=engine)

然后我就可以使用我的 x_rvs。问题是，我的发行版都不同。

一种选择是使用 [0,1] 范围内的正态分布并将其映射到其他分布。然而，考虑到 QMC 文档大力支持不同的发行版，我不禁觉得有一种方法可以直接生成我需要的随机变量。我知道我不能仅仅通过使用 QMC Sobol 引擎并重复使用它来生成不同的分布来生成准随机数据（如果我这样做，我得到的序列太相关并且没有用），所以我知道我必须在我的引擎中使用多个维度。

有人知道如何生成具有不同分布的多个准随机变量以保持 QMC 的良好特性吗？

Answer 1

但是，鉴于 QMC 文档非常重视支持不同的发行版......

此时，

scipy.stats.qmc

不具备从均匀分布、多项式和多元正态分布以外的任何分布直接生成 QMC 序列的功能。

你已经在

qrvs

中找到了类的

scipy.stats.sampling方法；这是 SciPy 在现有版本中唯一拥有的其他功能。

此外，SciPy 1.15 将拥有新的随机变量基础设施。您将能够使用

make_distribution

将（几乎）任何现有的 SciPy Continuous Distribution 转换为新接口，其中

sample

方法将接受 NumPy

Generator

或 SciPy

QMCEngine

作为

rng 的参数

参数。当它收到

QMCEngine

时，它会执行逆变换采样以从其他分布生成 QMC 样本。 （我相信这就是您所说的“一个选项是在 [0,1] 范围内使用[均匀]分布并将它们映射到其他分布”。）SciPy 1.15 RC1 现在可用。

我知道我不能仅仅通过使用 QMC Sobol 引擎并重复使用它来生成不同的分布来生成准随机数据（如果我这样做，我得到的序列太相关并且没有用）

是的，最好不要陷入那个陷阱！

所以我知道我必须在我的引擎中使用多个维度。

这是一个选项 - 您可以创建一个多元统一，在总和中每一项只有一个维度，然后使用每个维度的坐标作为“独立”一维、低差异样本。

另一种方法是从具有不同种子和

QMCEngine

的新

Generator

创建每个

scramble=True

。为此，SciPy 有一个私有

_rng_spawn

函数。它在

qmc_quad

中用于获取（对于大多数实际目的）积分的独立估计，并从中估计误差。 这里对这种方法进行了讨论。

所以在 SciPy 1.14 中可能看起来像：

import numpy as np
from scipy import stats

def _rng_spawn(rng, n_children):
    # https://github.com/scipy/scipy/blob/d916ddbab38468312fd0bee53ff68d5276f4f84f/scipy/_lib/_util.py#L1025
    # spawns independent RNGs from a parent RNG
    bg = rng._bit_generator
    ss = bg._seed_seq
    child_rngs = [np.random.Generator(type(bg)(child_ss))
                  for child_ss in ss.spawn(n_children)]
    return child_rngs

size = 16
rng = np.random.default_rng(235855178283348532492417763208635684641)  # high entropy seed generated with np.random.SeedSequence
rngs = _rng_spawn(rng, 2)

engine = stats.qmc.Sobol(d=1, scramble=True, seed=rngs[0])
dist = stats.fisk(3.9)
base = stats.sampling.NumericalInverseHermite(dist)
x_rvs1 = base.qrvs(size=size, qmc_engine=engine)

engine = stats.qmc.Sobol(d=1, scramble=True, seed=rngs[1])
dist = stats.norm()
base = stats.sampling.NumericalInverseHermite(dist)
x_rvs2 = base.qrvs(size=size, qmc_engine=engine)

x_rvs = x_rvs1 + x_rvs2