Was ist die effizienteste und portabel Möglichkeit, Gaußsche Zufallszahlen in Cython zu generieren?

Question

Ich schreibe eine Cython-Anwendung, wo ich eine Gaußsche Zufallsvariable on-the-fly in einer engen geschachtelten Schleife generieren muss. Ich möchte dies tun, ohne irgendwelche zusätzlichen Abhängigkeiten einzuführen, z. B. bei GSL.

Für eine Minimalversion der Art und Weise bin ich derzeit in der Lage dies mit einheitlichen Zufall Zahlen on-the-fly zu tun:

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 

Der obige Code funktional äquivalent zu numpy.random.rand ist (n) und kann mit der folgenden minimalen Setup-Datei zu erstellen:

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 

diese Frage zu beantworten, was ich suche ist die gleiche Art von Minimallösung für das funktionelle Äquivalent von np.random.randn (n), wieder ideal mit jedem Die Abhängigkeit wurde aus Gründen der Portabilität direkt aus libc.stdlib importiert.

Es gibt eine Beispielimplementierung auf the Wikipedia entry for the Box-Muller algorithm, aber ich hatte Schwierigkeiten, es wegen der Art zu implementieren, wie das konstante epsilon definiert wird.

Answer 1

Ich habe eine Funktion erstellt, die auf der Basis der polaren Version der Box-Muller-Transformation Gauß-verteilte Zufallszahlen erzeugt, wie im Pseudocode here beschrieben. (Die Wikipedia-Seite beschreibt diese Version ebenfalls, bietet aber keinen praktischen Pseudocode.)

Diese Methode generiert zwei Gauß-verteilte Zufallszahlen gleichzeitig. Das bedeutet, um die volle Geschwindigkeit zu erreichen, müssen wir einen Weg finden, um zwei Zahlen zu umgehen, ohne sie in Python-Objekte umzuwandeln. Der einfachste Weg dies zu tun (was ich mir vorstellen kann) ist, den Puffer zur direkten Manipulation durch den Generator zu übergeben. Das ist, was my_gaussian_fast tut, und es schlägt numpy mit einem bescheidenen Vorsprung.

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 

Tests. Hier ist eine einheitliche 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 

Das ist also auf jeden Fall schneller als numpy für normale Zufallszahl. Und wenn wir über sie schlau sind, dann ist es schneller für Gauß'sche Zufallszahl zu:

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 

Wie Robert Kern bestätigt, numpy beiden Werte verwendet erzeugt. wirft einen weg; my_gaussian_fast verwendet beide und speichert sie schnell.(Siehe die Geschichte dieser Antwort für eine naive my_gaussian_pair, die versucht, das Paar auf eine langsame Weise zurückzugeben.)

Answer 2

Sie sagen, dass Sie Probleme haben, die Umsetzung der Box-Muller wegen der Art und Weise transformieren sie definieren espilon:

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

According to here ist dies die C-Äquivalent:

const double lowest_double = -DBL_MAX; 

Also die richtige zu bekommen importiert in Cython:

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

Jetzt die Problem mit epsilon sollte gelöst werden.