Numpy гауссов из / dev / urandom - PullRequest
1 голос
/ 29 января 2020

У меня есть приложение, где мне нужно numpy.random.normal, но из источника crypgoraphi c PRNG. Numpy, по-видимому, не предоставляет эту опцию.

Лучшее, что я смог найти, было numpy.random.entropy.random_entropy, но это только uint32 и оно глючит, с большими массивами вы получаете "RuntimeError: Невозможно прочитать из системной криптографии c провайдер ", хотя urandom не является блокирующим ...

Однако вы можете сделать это: np.frombuffer(bytearray(os.urandom(1000*1000*4)), dtype=np.uint32).astype(np.double).reshape(1000, 1000)

Но у меня все еще остается проблема как-то преобразовать его в Гасианец и ничего не напортачил.

Есть ли решение, которое кто-то знает? Google отравлен numpy отбором из / dev / urandom, мне не нужно отсеивать, мне нужен urandom, являющийся единственным источником всей случайности.

1 Ответ

0 голосов
/ 29 января 2020

Я наткнулся на scipy.stats.rvs_ratio_uniforms и адаптировал их код для моих целей. Это всего в 3 раза медленнее, чем np.random.normal, несмотря на выборку в два раза больше случайности из источника c cryptographi.

import numpy as np
import os


def uniform_0_1(size):
    return np.frombuffer(bytearray(os.urandom(size*4)), dtype=np.uint32).astype(np.float) / 2**32

def normal(mu, sigma, size):
    bmp = np.sqrt(2.0/np.exp(1.0)) # about 0.8577638849607068

    size1d = tuple(np.atleast_1d(size))
    N = np.prod(size1d)  # number of rvs needed, reshape upon return

    x = np.zeros(N)
    simulated = 0

    while simulated < N:
        k = N - simulated

        a = uniform_0_1(size=k)
        b = (2.0 * uniform_0_1(size=k) - 1.0) * bmp

        accept = (b**2 <= - 4 * a**2 * np.log(a))
        num_accept = np.sum(accept)

        if num_accept > 0:
            x[simulated : (simulated + num_accept)] = (b[accept] * sigma / a[accept]) + mu
            simulated += num_accept


    return np.reshape(x, size1d)

Хотя одно беспокойство, numpy.random.random_sample: возвращать случайные числа с плавающей точкой в ​​полуоткрытом интервале [0.0 , 1.0).

Я не уверен, как добиться этой гарантии (никогда не 1,0) с моимiform_0_1 или даже имеет ли это значение.

...