Python/SciPy中关于低差异(例如Sobol)准随机序列的建议？

Question

我想在基于SciPy的模拟中使用准随机序列,特别是Sobol.关于现有高效套餐的任何建议？

Answer 1

是这个项目有什么用处？如果没有,相关的维基百科文章提到了NAG和Numerical Recipes中的C和Fortran例程,这可能不难在Python中包装或重新实现.GSL中还有C例程.

Answer 2

我会使用OpenTURNS，它提供了几个低差异序列：

• 福尔序列，
• 霍尔顿序列，
• 逆霍尔顿序列，
• 哈塞尔格罗夫序列，
• 索博尔序列。

此外，可以生成序列使得边缘具有任意分布。这是通过基于逆分布函数的概率变换来完成的。

在下面的示例中，我根据LowDiscrepancyExperiment类生成二维 Sobol' 序列。边际在 [-1, 1] 区间内均匀分布（这是 OT 中的默认均匀分布）。我建议使用等于 2 的幂的样本大小，因为 Sobol' 序列基于以 2 为基数的整数分解。该generate方法返回一个ot.Sample.

import openturns as ot
dim = 2
distribution = ot.ComposedDistribution([ot.Uniform()]*dim)
bounds = distribution.getRange()
sequence = ot.SobolSequence(dim)
samplesize = 2**5 # Sobol' sequences are in base 2
experiment = ot.LowDiscrepancyExperiment(sequence, distribution,
                                         samplesize, False)
sample = experiment.generate()
print(samplesize[:5])

前面的样本大小为 32。前 5 个元素是：

    y0      y1
0    0       0
1    0.5    -0.5
2   -0.5     0.5
3   -0.25   -0.25
4    0.75    0.75

OT 中的 Sobol' 序列可以生成任意大小的样本，维度高达 1111。

再做一点工作，我们就可以绘制出设计图。

import openturns.viewer as otv
fig = otv.PlotDesign(sample, bounds, 2**2, 2**1);
fig.set_size_inches(6, 6)

其产生：

看看每个基本区间中恰好有 4 个点。

如果需要，sample可以轻松转换为 Numpy 数组，这可能更适合您的 Scipy 要求：

import numpy as np
array = np.array(sample)

其他示例位于：http://openturns.github.io/openturns/master/examples/reliability_sensitivity/design_of_experiments.html