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 

はこの質問に答えるために、私は何を探していますと、n個np.random.randnと同等の機能のための最小限のソリューションと同じ種類の（あります）、もう一度理想的には移植性の理由から直接libc.stdlibに依存していました。

the Wikipedia entry for the Box-Muller algorithmの実装例がありますが、一定のイプシロンが定義されているため、実装に問題がありました。

あなたはので、彼らは espilonを定義する方法の変換ボックス・ミュラーを実装する問題があるとし

Answer 1

擬似コードhereで説明されているように、Box-Muller変換の極座標に基づいてガウス分布乱数を生成する関数を作成しました。 （ウィキペディアのページでもこのバージョンについて説明していますが、便利な擬似コードは提供していません）

2度にの乱数を生成します。つまり、完全なcythonの速度を得るには、2つの数値を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 * 2 

@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よりも高速です。我々はそれについてスマートなら、それはあまりにもガウス乱数のための高速です：Robert Kernによって確認されたように

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が1つ離れています。 my_gaussian_fastは両方を使用し、それらをすばやく保存します。（遅い方法でペアを返そうとしている純粋なmy_gaussian_pairについては、この回答の履歴を参照してください。）

Answer 2

：これはCと同等である

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

According to hereを：

const double lowest_double = -DBL_MAX; 

だから、正しいを取得しますCythonでのインポート：

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

epsilonの問題を解決する必要があります。