问题描述
我正在编写一个cython应用程序，需要在紧密的嵌套循环中即时生成一个高斯随机变量.我想这样做而不会引入任何额外的依赖关系，例如，对GSL的依赖.
I am writing a cython application where I need to generate a Gaussian random variable on-the-fly in a tight nested loop. I would like to do this without introducing any extra dependencies, e.g., on GSL.
对于当前我可以使用均匀随机数字即时实现的方法的最低版本:
For a minimal version of the way I am currently able to do this with uniform random numbers on-the-fly:
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)，并且可以使用以下最小安装文件进行编译:
The above code is functionally equivalent to numpy.random.rand(n), and can be compiled with the following minimal setup file:
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中导入任何依赖项.
To answer this question, what I am looking for is the same kind of minimal solution for the functional equivalent of np.random.randn(n), again ideally with any dependency directly cimported from libc.stdlib for reasons of portability.
Box-Muller算法的Wikipedia条目，但是由于定义了常数epsilon的方式，我在实现它时遇到了麻烦.
There is an example implementation on the Wikipedia entry for the Box-Muller algorithm, but I have had trouble implementing it due to the way that the constant epsilon is defined.
推荐答案
我创建了一个函数，该函数根据Box-Muller变换的极坐标版本生成高斯分布的随机数，如伪代码此处.(我最初在存档此处.)
I created a function that generates gaussian-distributed random numbers based on the polar version of the Box-Muller transformation, as described by the pseudocode here. (I originally found this on the page archived here.)
此方法一次生成两个高斯分布的随机数.这意味着要获得完整的 cython 速度，我们需要找到一种方法来传递两个数字而不将其转换为Python对象.(我能想到的)最简单的方法是将缓冲区传递给生成器直接操作.这就是 my_gaussian_fast 所做的，并且以适度的优势击败了 numpy .
This method generates two gaussian-distributed random numbers at a time. That means to get full cython speed, we need to figure out a way to pass two numbers around without turning them into Python objects. The most straightforward way to do so (that I can think of) is to pass the buffer in for direct manipulation by the generator. That's what my_gaussian_fast does, and it beats numpy by a modest margin.
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
测试.这是一个统一的基准:
Tests. Here's a uniform benchmark:
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 快.而且，如果我们对此很聪明，那么对高斯随机数也更快:
So this is definitely faster than numpy for ordinary random numbers. And if we're smart about it, it's faster for gaussian random numbers too:
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 同时使用和存储它们.(请参阅此答案的历史记录，以尝试以慢速方式返回对的天真 my_gaussian_pair .)
As confirmed by Robert Kern, numpy uses both values generated. my_gaussian throws one away; my_gaussian_fast uses both and stores them quickly. (See this answer's history for a naive my_gaussian_pair that tries to return the pair in a slow way.)
这篇关于在cython中生成高斯随机数的最有效，最便捷的方法是什么?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！