Как эффективно выполнять миллиарды извлечений Бернулли с помощью numpy?

Question

Я работаю над диссертацией по эпидемиологии, и мне нужно смоделировать эпидемию SI во временной сети. На каждом временном шаге есть вероятность ~ Бернулли (бета) выполнить извлечение между зараженным и уязвимым узлом. Я использую np.random.binomial(size=any, n=1, p=beta), чтобы компьютер решил. Теперь мне нужно смоделировать эпидемию в той же сети, запустив ее с каждого из узлов. Это следует повторить K раз, чтобы получить некоторые статистически значимые результаты для каждого узла, и, поскольку временная сеть также является стохастической, все должно быть повторено NET_REALIZATION раз.

Итак, в сети с N = 100, если K=500 и NET=REALIZATION=500, эпидемия должна повториться 25 000 000 раз. Если T=100, это означает 2500 000 000 извлечений на набор пар SI (что, конечно, меняется во времени). Если бета-версия мала, что часто бывает, это приводит к очень трудоемким вычислениям. Если вы думаете, что для моего компьютера извлечение Бернулли занимает 3,63 мкс, это означает, что мне нужно ждать несколько часов, чтобы получить какие-то результаты, что действительно ограничивает развитие моей диссертации. Проблема в том, что больше половины времени тратится на случайные извлечения. Я должен использовать numpy, поскольку результаты извлечения взаимодействуют с другими структурами данных. Я пробовал использовать numba, но, похоже, это не улучшило скорость извлечения.Есть ли более быстрый способ получить такие же результаты? Я думал о том, чтобы сделать очень-очень большое извлечение один раз навсегда, что-то вроде извлечения 10^12 нулей и единиц, и просто импортировать их часть для каждой другой симуляции (это следует повторить для нескольких значений бета), но мне интересно если есть более умный ход.

Спасибо за помощь

Answer 1

Если вы можете выразить свое betas как приращение 2^-N (например, приращение 1/256, если N равно 8.), затем извлеките случайные 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 представляет новую систему генерации случайных чисел наряду с устаревшей. Возможно, он имеет лучшую производительность, генерирующую Бернулли и биномиальные переменные, чем старый, особенно потому, что его ГСЧ по умолчанию, PCG64, легче, чем Мерсенн Твистер по умолчанию в устаревшей системе. Ниже приводится пример.

import numpy

beta = 5.0/256
rs = numpy.random.default_rng()
rbinom = rs.binomial(10, beta)