在 cython 中生成高斯随机数的最有效和可移植的方法是什么？

Question

我正在编写一个 cython 应用程序，我需要在一个紧密的嵌套循环中即时生成一个高斯随机变量。我想在不引入任何额外依赖项的情况下执行此操作，例如，在 GSL 上。

对于我目前能够使用 均匀随机 动态数字执行此操作的最小版本：

from libc.stdlib cimport rand, RAND_MAX
import numpy as np

cdef double random_uniform():
    cdef double r = rand()
    return r/RAND_MAX

def my_function(int n):
    cdef int i
    cdef double[:] result = np.zeros(n, dtype='f8', order='C')
    for i in range(n):
        result[i] = random_uniform()
    return result

以上代码在功能上等同于numpy.random.rand(n)，可以使用以下最小安装文件进行编译：

from distutils.core import setup
from Cython.Build import cythonize
import numpy as np

setup(ext_modules=cythonize("example.pyx"), include_dirs=[np.get_include()])

# compile instructions:
# python setup.py build_ext --inplace

为了回答这个问题，我正在寻找的是与 np.random.randn(n) 的功能等效的同一种最小解决方案，同样理想情况下，任何依赖直接从 libc.stdlib 导入便携性的原因。

the Wikipedia entry for the Box-Muller algorithm 上有一个示例实现，但由于常量 epsilon 的定义方式，我在实现它时遇到了麻烦。

Answer 1

你说你在实施 Box-Muller 变换时遇到问题，因为它们定义的方式 espilon:

const double epsilon = std::numeric_limits<double>::min();

According to here 这是 C 等价物：

const double lowest_double = -DBL_MAX;

所以要在 Cython 中获得正确的导入：

from libc.float import DBL_MAX #it should still be portable btw.

现在 epsilon 的问题应该已经解决了。

Answer 2

我创建了一个函数，它根据 Box-Muller 变换的极坐标版本生成高斯分布的随机数，如伪代码 here. (I originally found this on the page archived here 所述。)

此方法一次生成两个个高斯分布的随机数。这意味着要获得完整的 cython 速度，我们需要找出一种方法来传递两个数字而不将它们变成 Python 对象。最直接的方法（我能想到的）是将缓冲区传递给生成器直接操作。这就是 my_gaussian_fast 所做的，它以适度的优势击败了 numpy。

from libc.stdlib cimport rand, RAND_MAX
from libc.math cimport log, sqrt
import numpy as np
import cython

cdef double random_uniform():
    cdef double r = rand()
    return r / RAND_MAX

cdef double random_gaussian():
    cdef double x1, x2, w

    w = 2.0
    while (w >= 1.0):
        x1 = 2.0 * random_uniform() - 1.0
        x2 = 2.0 * random_uniform() - 1.0
        w = x1 * x1 + x2 * x2

    w = ((-2.0 * log(w)) / w) ** 0.5
    return x1 * w

@cython.boundscheck(False)
cdef void assign_random_gaussian_pair(double[:] out, int assign_ix):
    cdef double x1, x2, w

    w = 2.0
    while (w >= 1.0):
        x1 = 2.0 * random_uniform() - 1.0
        x2 = 2.0 * random_uniform() - 1.0
        w = x1 * x1 + x2 * x2

    w = sqrt((-2.0 * log(w)) / w)
    out[assign_ix] = x1 * w
    out[assign_ix + 1] = x2 * w

@cython.boundscheck(False)
def my_uniform(int n):
    cdef int i
    cdef double[:] result = np.zeros(n, dtype='f8', order='C')
    for i in range(n):
        result[i] = random_uniform()
    return result

@cython.boundscheck(False)
def my_gaussian(int n):
    cdef int i
    cdef double[:] result = np.zeros(n, dtype='f8', order='C')
    for i in range(n):
        result[i] = random_gaussian()
    return result

@cython.boundscheck(False)
def my_gaussian_fast(int n):
    cdef int i
    cdef double[:] result = np.zeros(n, dtype='f8', order='C')
    for i in range(n // 2):  # Int division ensures trailing index if n is odd.
        assign_random_gaussian_pair(result, i * 2)
    if n % 2 == 1:
        result[n - 1] = random_gaussian()

    return result

测试。这是一个统一的基准：

In [3]: %timeit numpy.random.uniform(size=10000)
10000 loops, best of 3: 130 µs per loop

In [4]: %timeit numpy.array(example.my_uniform(10000))
10000 loops, best of 3: 85.4 µs per loop

所以对于普通的随机数来说，这肯定比numpy快。如果我们对此很聪明，高斯随机数也会更快：

In [5]: %timeit numpy.random.normal(size=10000)
1000 loops, best of 3: 393 µs per loop

In [6]: %timeit numpy.array(example.my_gaussian(10000))
1000 loops, best of 3: 542 µs per loop

In [7]: %timeit numpy.array(example.my_gaussian_fast(10000))
1000 loops, best of 3: 266 µs per loop

正如确认的那样，numpy 使用生成的两个值。 my_gaussian 扔掉一个； my_gaussian_fast 同时使用并快速存储它们。（请参阅此答案的历史，以了解试图以缓慢的方式 return 这对天真的 my_gaussian_pair。）

在 cython 中生成高斯随机数的最有效和可移植的方法是什么？

What is the most efficient and portable way to generate Gaussian random numbers in cython?

python

performance

numpy

cython

scientific-computing