来自 /dev/urandom 的 Numpy Gaussian答案

【问题标题】：Numpy Gaussian from /dev/urandom来自 /dev/urandom 的 Numpy Gaussian
【发布时间】：2020-05-15 05:24:27
【问题描述】：

我有一个需要 numpy.random.normal 的应用程序，但来自一个神秘的 PRNG 源。 Numpy 似乎没有提供这个选项。

我能找到的最好的是numpy.random.entropy.random_entropy，但这只是 uint32 并且有问题，即使 urandom 是非阻塞的，对于大型数组，您也会收到“RuntimeError: Unable to read from system cryptographic provider”...

您可以这样做：np.frombuffer(bytearray(os.urandom(1000*1000*4)), dtype=np.uint32).astype(np.double).reshape(1000, 1000)

但我仍然存在以某种方式将其转换为高斯而不是搞砸任何事情的问题。

有没有人知道的解决方案？ Google 被 /dev/urandom 的 numpy 播种毒害了，我不需要播种，我需要 urandom 作为所有随机性的唯一来源。

【问题讨论】：

我想通过适当的转换等，您可以从 dev/urandom 获得 uniform 分布。因此，您要查找的内容类似于“如何将均匀分布转换为正态分布”：mathoverflow.net/questions/28557/…
将从 urandom 采样的 int 值归一化为 [0,1]，然后使用逆 CDF 转换为正态分布的随机数。（更新：Lagerbaer 的类似建议相同）
为了补充@Lagerbaer 链接的答案，this answer 有一个极地 Box-Muller 变换的快速 cython 实现，您可能会发现它很有用。

标签： python numpy random gaussian normal-distribution

【解决方案1】：

我遇到了scipy.stats.rvs_ratio_uniforms，并根据我的目的修改了他们的代码。它只比np.random.normal 慢 3 倍，尽管从加密源中采样两倍的随机性。

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) 中返回随机浮点数。

我不确定如何使用我的 uniform_0_1 来实现该保证（从不 1.0），或者它是否重要。

【讨论】：