在 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 导入任何依赖项。

维基百科条目上有Box-Muller 算法的示例实现，但由于常数 epsilon 的定义方式，我在实现它时遇到了麻烦。

Answer 1

我创建了一个函数，该函数根据 Box-Muller 变换的极坐标版本生成高斯分布随机数，如此处伪代码所述。（我最初是在此处存档的页面上找到的。）

\n

该方法一次生成两个高斯分布的随机数。这意味着要获得全速cython，我们需要找到一种方法来传递两个数字而不将它们转换为 Python 对象。最直接的方法（我能想到的）是传递缓冲区以供生成器直接操作。事实就是如此my_gaussian_fast，而且它numpy以微弱优势击败了对手。

\n

from libc.stdlib cimport rand, RAND_MAX\nfrom libc.math cimport log, sqrt\nimport numpy as np\nimport cython\n\ncdef double random_uniform():\n    cdef double r = rand()\n    return r / RAND_MAX\n\ncdef double random_gaussian():\n    cdef double x1, x2, w\n\n    w = 2.0\n    while (w >= 1.0):\n        x1 = 2.0 * random_uniform() - 1.0\n        x2 = 2.0 * random_uniform() - 1.0\n        w = x1 * x1 + x2 * x2\n\n    w = ((-2.0 * log(w)) / w) ** 0.5\n    return x1 * w\n\n@cython.boundscheck(False)\ncdef void assign_random_gaussian_pair(double[:] out, int assign_ix):\n    cdef double x1, x2, w\n\n    w = 2.0\n    while (w >= 1.0):\n        x1 = 2.0 * random_uniform() - 1.0\n        x2 = 2.0 * random_uniform() - 1.0\n        w = x1 * x1 + x2 * x2\n\n    w = sqrt((-2.0 * log(w)) / w)\n    out[assign_ix] = x1 * w\n    out[assign_ix + 1] = x2 * w\n\n@cython.boundscheck(False)\ndef my_uniform(int n):\n    cdef int i\n    cdef double[:] result = np.zeros(n, dtype=\'f8\', order=\'C\')\n    for i in range(n):\n        result[i] = random_uniform()\n    return result\n\n@cython.boundscheck(False)\ndef my_gaussian(int n):\n    cdef int i\n    cdef double[:] result = np.zeros(n, dtype=\'f8\', order=\'C\')\n    for i in range(n):\n        result[i] = random_gaussian()\n    return result\n\n@cython.boundscheck(False)\ndef my_gaussian_fast(int n):\n    cdef int i\n    cdef double[:] result = np.zeros(n, dtype=\'f8\', order=\'C\')\n    for i in range(n // 2):  # Int division ensures trailing index if n is odd.\n        assign_random_gaussian_pair(result, i * 2)\n    if n % 2 == 1:\n        result[n - 1] = random_gaussian()\n\n    return result\n

\n

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

\n

In [3]: %timeit numpy.random.uniform(size=10000)\n10000 loops, best of 3: 130 \xc2\xb5s per loop\n\nIn [4]: %timeit numpy.array(example.my_uniform(10000))\n10000 loops, best of 3: 85.4 \xc2\xb5s per loop\n

\n

所以这肯定比numpy普通随机数要快。如果我们聪明的话，高斯随机数也会更快：

\n

In [5]: %timeit numpy.random.normal(size=10000)\n1000 loops, best of 3: 393 \xc2\xb5s per loop\n\nIn [6]: %timeit numpy.array(example.my_gaussian(10000))\n1000 loops, best of 3: 542 \xc2\xb5s per loop\n\nIn [7]: %timeit numpy.array(example.my_gaussian_fast(10000))\n1000 loops, best of 3: 266 \xc2\xb5s per loop\n

\n

正如Robert Kern所确认的，numpy使用生成的两个值。my_gaussian扔掉一个；my_gaussian_fast两者都可以使用并可以快速存储。（请参阅此答案的历史记录，了解my_gaussian_pair尝试以缓慢的方式返回该对的天真想法。）

\n