0

我正在尝试使用 python 进行准随机 Halton 采样和排序。有没有我可以使用的关于这个问题并绘制输出的库？如果您能帮助我，我将不胜感激。

4

2 回答 2

1

在 Python 中，SciPy 是主要的科学计算包，它包含一个 Halton 序列生成器，以及其他 QMC 函数。

对于绘图，SciPy 的标准方法是 matplotlib；如果您对此不熟悉，SciPy 教程也是一个很好的起点。

一个基本的例子：

from scipy.stats import qmc
import matplotlib.pyplot as plt

sampler = qmc.Halton(d=2, scramble=True)
sample = sampler.random(n=5)

plt.plot(sample)
于 2021-08-27T18:46:51.477 回答
0

如果您想使用第一原则来做到这一点，您可以使用Wikipedia 页面上提供的 Halton 序列代码：

def halton(b):
    """Generator function for Halton sequence."""
    n, d = 0, 1
    while True:
        x = d - n
        if x == 1:
            n = 1
            d *= b
        else:
            y = d // b
            while x <= y:
                y //= b
            n = (b + 1) * y - x
        yield n / d

这是一种获取 2 和 3 序列的前 256 个点并绘制它们的方法：

n = 256
df = pd.DataFrame([
    (x, y) for _, x, y in zip(range(n), halton(2), halton(3))
], columns=list('xy'))

ax = df.plot.scatter(x='x', y='y')
ax.set_aspect('equal')

附录“如果我想为所有序列取 200 个点和 8 个维度，我应该怎么做？”

假设您可以访问素数序列生成器，那么您可以使用例如：

m = 8
p = list(primes(m * m))[:m]  # first m primes

n = 200
z = pd.DataFrame([a for _, *a in zip(range(n), *[halton(k) for k in p])])

>>> z
            0         1       2         3         4         5         6         7
0    0.500000  0.333333  0.2000  0.142857  0.090909  0.076923  0.058824  0.052632
1    0.250000  0.666667  0.4000  0.285714  0.181818  0.153846  0.117647  0.105263
2    0.750000  0.111111  0.6000  0.428571  0.272727  0.230769  0.176471  0.157895
3    0.125000  0.444444  0.8000  0.571429  0.363636  0.307692  0.235294  0.210526
4    0.625000  0.777778  0.0400  0.714286  0.454545  0.384615  0.294118  0.263158
..        ...       ...     ...       ...       ...       ...       ...       ...
195  0.136719  0.576132  0.3776  0.011662  0.868520  0.089213  0.567474  0.343490
196  0.636719  0.909465  0.5776  0.154519  0.959429  0.166136  0.626298  0.396122
197  0.386719  0.057613  0.7776  0.297376  0.058603  0.243059  0.685121  0.448753
198  0.886719  0.390947  0.9776  0.440233  0.149512  0.319982  0.743945  0.501385
199  0.074219  0.724280  0.0256  0.583090  0.240421  0.396905  0.802768  0.554017

奖金

使用Atkin 筛的素数序列：

import numpy as np

def primes(limit):
    # Generates prime numbers between 2 and n
    # Atkin's sieve -- see http://en.wikipedia.org/wiki/Prime_number
    sqrtLimit = int(np.sqrt(limit)) + 1

    # initialize the sieve
    is_prime = [False, False, True, True, False] + [False for _ in range(5, limit + 1)]

    # put in candidate primes:
    # integers which have an odd number of
    # representations by certain quadratic forms
    for x in range(1, sqrtLimit):
        x2 = x * x
        for y in range(1, sqrtLimit):
            y2 = y*y
            n = 4 * x2 + y2
            if n <= limit and (n % 12 == 1 or n % 12 == 5): is_prime[n] ^= True
            n = 3 * x2 + y2
            if n <= limit and (n % 12 == 7): is_prime[n] ^= True
            n = 3*x2-y2
            if n <= limit and x > y and n % 12 == 11: is_prime[n] ^= True

    # eliminate composites by sieving
    for n in range(5, sqrtLimit):
        if is_prime[n]:
            sqN = n**2
            # n is prime, omit multiples of its square; this is sufficient because
            # composites which managed to get on the list cannot be square-free
            for i in range(1, int(limit/sqN) + 1):
                k = i * sqN # k ∈ {n², 2n², 3n², ..., limit}
                is_prime[k] = False
    for i, truth in enumerate(is_prime):
        if truth: yield i
于 2021-08-27T19:00:51.617 回答