如何使用 numpy 高效地执行数十亿次伯努利提取？
How to efficiently perform billions of Bernoulli extractions using numpy?
我正在写一篇关于流行病学的论文，我必须在时间网络中模拟 SI 流行病。在每个时间步都有概率 ~ Bernoulli(beta) 在受感染节点和易感节点之间执行提取。我正在使用 np.random.binomial(size=whatever, n=1, p=beta) 让计算机决定。现在，我必须通过从每个节点开始来模拟同一网络中的流行病。这应该重复 K 次以获得每个节点的一些统计相关结果，并且由于时间网络也是随机的，所以一切都应该重复 NET_REALIZATION 次。
所以，在一个N=100的网络中，如果K=500，NET=REALIZATION=500，流行病应该重复25,000,000次。如果 T=100，则意味着每组 S-I 对提取 2,500,000,000 次（当然会随时间变化）。如果 beta 很小（通常是这种情况），这会导致非常耗时的计算。
如果你认为，对于我的计算机，伯努利提取需要 3.63 微秒才能发生，这意味着我必须等待数小时才能得到一些结果，这确实限制了我论文的发展。
问题是一半以上的时间都花在了随机提取上。
我应该使用 numpy，因为提取的结果会与其他数据结构交互。我尝试使用 numba，但它似乎并没有提高提取速度。
有没有更快的方法来获得相同的结果？我正在考虑永远做一次非常非常大的提取，比如 0 和 1 的 10^12 提取，然后为每个不同的模拟导入其中的一部分（这应该针对多个 beta 值重复），但我想知道如果有更聪明的举动。
感谢帮助
如果您可以将 betas 表示为 2^-N 的增量（例如，如果 N 为 8，则增量为 1/256。），然后提取随机的 N 位块并确定每个块是否块小于 beta * 2^N。如果 32 能被 N 整除，这会更好。
请注意，numpy.random.uniform 生成随机浮点数，预计比生成随机整数或位慢。这尤其是因为生成随机浮点数取决于生成随机整数——而不是相反。
以下是这个想法如何运作的一个例子。
import numpy
# Fixed seed for demonstration purposes
rs = numpy.random.RandomState(777778)
# Generate 10 integers in [0, 256)
ri = rs.randint(0, 256, 10)
# Now each integer x can be expressed, say, as a Bernoulli(5/256)
# variable which is 0 if x < 5, and 1 otherwise. I haven't tested
# the following, which is similar to an example you gave in a
# comment.
rbern = (ri>=5) * 1
如果您可以使用 NumPy 1.17 或更高版本，则存在以下替代方案：
import numpy
rs = numpy.random.default_rng()
ri = rs.integers(0, 256, 10)
另请注意，NumPy 1.17 引入了一个 new random number generation system 以及旧版本。也许它在生成伯努利变量和二项式变量方面比旧版本具有更好的性能，尤其是因为它的默认 RNG PCG64 比遗留系统的默认 Mersenne Twister 更轻量级。下面是一个例子。
import numpy
beta = 5.0/256
rs = numpy.random.default_rng()
rbinom = rs.binomial(10, beta)
我正在写一篇关于流行病学的论文，我必须在时间网络中模拟 SI 流行病。在每个时间步都有概率 ~ Bernoulli(beta) 在受感染节点和易感节点之间执行提取。我正在使用 np.random.binomial(size=whatever, n=1, p=beta) 让计算机决定。现在，我必须通过从每个节点开始来模拟同一网络中的流行病。这应该重复 K 次以获得每个节点的一些统计相关结果，并且由于时间网络也是随机的，所以一切都应该重复 NET_REALIZATION 次。
所以，在一个N=100的网络中，如果K=500，NET=REALIZATION=500，流行病应该重复25,000,000次。如果 T=100，则意味着每组 S-I 对提取 2,500,000,000 次（当然会随时间变化）。如果 beta 很小（通常是这种情况），这会导致非常耗时的计算。 如果你认为，对于我的计算机，伯努利提取需要 3.63 微秒才能发生，这意味着我必须等待数小时才能得到一些结果，这确实限制了我论文的发展。 问题是一半以上的时间都花在了随机提取上。 我应该使用 numpy，因为提取的结果会与其他数据结构交互。我尝试使用 numba，但它似乎并没有提高提取速度。 有没有更快的方法来获得相同的结果？我正在考虑永远做一次非常非常大的提取，比如 0 和 1 的 10^12 提取，然后为每个不同的模拟导入其中的一部分（这应该针对多个 beta 值重复），但我想知道如果有更聪明的举动。
感谢帮助
如果您可以将 betas 表示为 2^-N 的增量（例如，如果 N 为 8，则增量为 1/256。），然后提取随机的 N 位块并确定每个块是否块小于 beta * 2^N。如果 32 能被 N 整除，这会更好。
请注意，numpy.random.uniform 生成随机浮点数，预计比生成随机整数或位慢。这尤其是因为生成随机浮点数取决于生成随机整数——而不是相反。
以下是这个想法如何运作的一个例子。
import numpy
# Fixed seed for demonstration purposes
rs = numpy.random.RandomState(777778)
# Generate 10 integers in [0, 256)
ri = rs.randint(0, 256, 10)
# Now each integer x can be expressed, say, as a Bernoulli(5/256)
# variable which is 0 if x < 5, and 1 otherwise. I haven't tested
# the following, which is similar to an example you gave in a
# comment.
rbern = (ri>=5) * 1
如果您可以使用 NumPy 1.17 或更高版本，则存在以下替代方案：
import numpy
rs = numpy.random.default_rng()
ri = rs.integers(0, 256, 10)
另请注意，NumPy 1.17 引入了一个 new random number generation system 以及旧版本。也许它在生成伯努利变量和二项式变量方面比旧版本具有更好的性能，尤其是因为它的默认 RNG PCG64 比遗留系统的默认 Mersenne Twister 更轻量级。下面是一个例子。
import numpy
beta = 5.0/256
rs = numpy.random.default_rng()
rbinom = rs.binomial(10, beta)