1

我正在写一篇关于流行病学的论文，我必须在时间网络中模拟 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，则意味着每组 SI 对有 2,500,000,000 次提取（当然随时间变化）。如果 beta 很小（通常是这种情况），这会导致计算非常耗时。如果你认为，对于我的电脑来说，伯努利提取需要 3.63 µs，这意味着我必须等待几个小时才能得到一些结果，这确实限制了我的论文的发展。问题是超过一半的时间都花在了随机抽取上。我应该使用 numpy，因为提取的结果与其他数据结构交互。我尝试使用 numba，但它似乎并没有提高提取速度。有没有更快的方法来获得相同的结果？我正在考虑永远做一次非常大的提取，比如 0 和 1 的 10 ^ 12 次提取，然后为每个不同的模拟导入其中的一部分（这应该针对几个 beta 值重复），但我想知道如果有更聪明的举动。

感谢帮助

4

1 回答 1

1

如果您可以将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 引入了一个新的随机数生成系统以及旧系统。也许它在生成伯努利和二项式变量时比旧的性能更好，特别是因为它的默认 RNG PCG64 比旧系统的默认 Mersenne Twister 更轻。下面是一个例子。

import numpy

beta = 5.0/256
rs = numpy.random.default_rng()
rbinom = rs.binomial(10, beta)
于 2019-11-24T10:32:56.460 回答