logo标签列表

在Python/Numpy/Scipy中生成低偏差准随机序列？

pythonnumpyrandomnumbersscipy
19

这个问题已经有一个了，但是答案包含了一个错误的链接，而且已经超过两年了，我希望现在有更好的解决方案 :)

低差异准随机序列，例如Sobol序列，比均匀随机序列更加均匀地填充空间。是否有一种好的/简单的方法在Python中生成它们？

4个回答
17

我认为在Python中替代低差异序列的最佳选择是敏感性分析库（SALib）：

https://github.com/SALib/SALib

我认为这是一个活跃的项目，您可以联系作者检查您需要的功能是否已经实现。如果这不能解决您的问题，Corrado Chisari 将John Burkardt在Matlab中制作的SOBOL版本移植到了Python中，您可以在此处访问：

http://people.sc.fsu.edu/~jburkardt/py_src/sobol/sobol.html

有人清理了这些源代码中的注释，并将它们放入了文档字符串的格式中。它更易读，您可以在此处访问：

https://github.com/naught101/sobol_seq

您列出的sobol-seq软件包将维度限制在40。您是否知道有哪些软件包可以生成更高维度的Sobol数？ - toing
4
scipy.stats.qmc 在最新的 scipy 发行版（1.6.2）中似乎不可用。 - G. Fougeron
2
我认为现在（SciPy版本>=1.7.1），最简单的方法就是我在这里所做的方式。他们实施了Joe和Kuo算法，可处理高达21,201个维度，这是你可以获得的最高维数（开源）。https://web.maths.unsw.edu.au/~fkuo/sobol/ 在这里，我展示了如何使用base2方法（带有Owen混淆）和随机方法（从序列生成任意数量的点），以及如何跳过第一个点。
请注意，这个例程可能会非常慢（由于将点转换为冲击的ndtri或反向正态分布转换），特别是在高维度+高模拟计数的情况下。从Sobol序列生成点非常快，但对于大多数蒙特卡罗模拟，您需要将它们转换为冲击（您可能正在使用标准正态分布之外的其他分布）。这至少让您能够直接在Python代码中生成点。
此外，在QMCgenerate例程中，我跳过了第一个点（即0秒）-虽然这是常见的做法，但一些论文建议不要这样做（但如果您有更好的替代方案，请随时评论）。我之所以转置它们，只是为了稍后将其粘贴到Excel中并检查生成的冲击。希望那些需要这个算法的人会觉得有用。
from scipy.stats import qmc # needs SciPy >= 1.7.1
from scipy.special import ndtri
import numpy as np
import timeit

time_periods = 252
factors = 12

# IF using base2 generation, need a pow(2,m)
sims = 8192 

dimensions = factors*time_periods

def RQMCgenerate (dimensions, sims, seed):
    start_time = timeit.default_timer()
    m=10 # start at 1024 sims
    while pow(2,m) < sims: #m = 17 # 131,072 sims; M = 16 # 65,536 sims
        m = m+1
    RQMCgenerator = qmc.Sobol(dimensions, scramble=True, seed=seed)
    RQMCsamples = RQMCgenerator.random_base2(m)
    print('\n' + 'Time after sample generation RQMC:', (timeit.default_timer() - start_time), 'seconds'); 
    sobol = ndtri(RQMCsamples).T # get normsinv(points) and transpose to dimensions * sims 
    del RQMCsamples
    print('\n' + 'Time after ndtri (normsinv) of', sims,'sims x dimensions', dimensions, 'Randomized Sobol points): ', (timeit.default_timer() - start_time), 'seconds');
    return sobol

def QMCgenerate(dimensions, sims):
    start_time = timeit.default_timer()
    QMCgenerator = qmc.Sobol(dimensions, scramble=False)
    QMCgenerator.fast_forward(1) #skip first point where normsinv(0) = -Inf
    QMCsamples = QMCgenerator.random(sims) #this generates points not having to be powers of 2
    print('\n' + 'Time after sample generation QMC:', (timeit.default_timer() - start_time), 'seconds'); 
    sobol = ndtri(QMCsamples).T # get normsinv(points) and transpose to dimensions * sims
    del QMCsamples
    print('\n' + 'Time after ndtri (normsinv) of', sims,'sims x dimensions', dimensions, 'Sobol points):', (timeit.default_timer() - start_time), 'seconds');
    return sobol

RQMCsobol = RQMCgenerate(dimensions, sims, seed=0) #note sims changed with pow(2,m) if a power of 2 was not passed
sobol = QMCgenerate(dimensions, sims)

样本生成后的RQMC时间：0.4269224999952712秒

8092个模拟x 3024个随机化Sobol点的ndtri（normsinv）之后的时间：1.0048970999996527秒

样本生成后的QMC时间：0.0630135999963386秒

8092个模拟x 3024个Sobol点的ndtri（normsinv）之后的时间：0.5444753999981913秒

当sims * dimensions变得更高时，速度会变得更慢，尽管我还没有在Python中找到比ndtri更快的将点转换为正态分布冲击的方法：

样本生成后的RQMC时间：2.1779929000040283秒

131072个模拟x 3024个随机化Sobol点的ndtri（normsinv）之后的时间：10.617904700004146秒

样本生成后的QMC时间：1.079756200000702秒

131072个模拟x 3024个Sobol点的ndtri（normsinv）之后的时间：9.545934699999634秒

SciPy的维护者和QMC子模块的作者在此。请不要跳过第一点。我们与所有QMC社区进行了广泛的讨论，发现这样做只会带来缺点。如果您不想有一些0，请使用混淆。Art Owen撰写了一篇关于此的论文https://arxiv.org/abs/2008.08051。 - tupui
对于使用非均匀分布的抽样，最好使用反变换（例如请参见MultivariateNormalQMC）或者使用我们提供的处理QMC随机变量的抽样方法（即将在1.9版本中推出）。http://scipy.github.io/devdocs/reference/generated/scipy.stats.sampling.NumericalInverseHermite.html - tupui
@tupui 我意识到关于QMC序列的研究论文已经确定跳过第一个点是不推荐的；然而，在实践中，仍有许多人跳过它。如果你有足够的模拟，那么从实用的角度来看，这并不重要。 - Matt
没有关系，这正是文章的重点。请看数据。只要缺少一个点，性质就会丧失，收敛速度会急剧下降。而且随着点数的增加，也无法恢复。Sobol'是2的幂序列，只有在特定条件下才能像广告中所说的那样运行。否则，您只会得到普通的QMC速率1/n，而不是1/(3/2)。 - tupui
本文展示了如何使用跳跃参数以及生成的点分布的变化。在期权定价中跳过一个点几乎没有影响，除非是学术界认为有影响。第一个点是0，这意味着如果您从反向正态分布或任何反向分布中提取数据，那么就不会有任何变化。如果使用的点数为1024，则最大误差加权为1/1024，约为0.000977。我只是在谈论现实，而不是关于某些数学上不正确但在实践中很好的学术练习。 - Matt
显示剩余2条评论
1

Chaospy也是一种有效的选择。可以选择多种低差异采样方法（包括'Sobol'、'latin hypercube'等）- 更多详情请参见文档。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接