Эффективное моделирование множества частотных распределений за тысячи итераций в numpy ⇐ Python

[Anonymous]

У меня на работе возникла следующая проблема:
У нас есть, скажем, 1 миллион возможных событий, которые определяют распределение частоты и серьезности.
Для каждого события у нас есть годовая ставка, которая определяет распределение Пуассона, а также параметры альфа и бета для бета-распределений. Цель состоит в том, чтобы смоделировать порядка >100 000 «лет», при этом каждый год определяется как получение частоты N для каждого события и получение N выборок относительного бета-распределения.
Обременительным для меня фактом является то, как я могу эффективно получить образцы N_i ~ Пуассона (lambda_i) из бета-распределения Beta_i, одновременно гарантируя, что я могу отнести их к правильному году?
Что касается результатов, я мне нужно будет посмотреть как максимальное, так и общее значение образцов за год, поэтому временно я просто сохраняю его
как массив словарей (не предназначенный для использования в качестве выходного формата)

Код: Выделить всё

years = 5000

rng = np.random.default_rng()
losses = []
for year in range(years):
occurences = rng.poisson(data['RATE'])
annual_losses = []
for idx, occs in enumerate(occurences):
if occs > 0:
event = data.iloc[idx]
for occ in range(occs):
loss = rng.beta(event['alpha'], event['beta']) * event['ExpValue']
annual_losses.append(loss)
annual_losses.append(0)
losses.append({'year': year, 'losses': annual_losses})

Я пытался выполнить оптимизацию кода Python/Numpy, используемого для моделирования, но не могу понять, как эффективно векторизовать этот код.
Изменения, которые я сделал перед публикацией здесь (раз за 5000 лет):

• замена с scipy на numpy (72 с -> 66 с)
• вычисление частот за все годы за один раз вне цикла (66 с -> 73 с... упс)

В идеале я бы вот так, чтобы работать как можно быстрее или выполнять как можно больше итераций, и ранее у меня также возникали проблемы с памятью при использовании scipy.
РЕДАКТИРОВАТЬ:
По запросу вот версия, в которой просто рассчитывается максимальный убыток за каждый год

Код: Выделить всё

years = 20_000

rng = np.random.default_rng()
largest_losses = np.zeros(shape=years)
for year in range(years):
occurences = rng.poisson(data['RATE'])
largest_loss = 0
for idx, event_occs in enumerate(occurences):
if event_occs > 0:
event = data.iloc[idx]
for event_occ in range(event_occs):
loss = rng.beta(event['alpha'], event['beta']) * event['ExpValue']
if loss > largest_loss:
largest_loss = loss
largest_losses[year] = largest_loss

Для целей тестирования общая частота событий составляет ~0,997, а приведенные выше тесты по времени взяты из пула событий, состоящего из 100 604 событий.
Чтобы указать цель, я хотел бы посмотреть, смогу ли я рассчитать процентили потерь, например. потеря 1 из 250

Код: Выделить всё

np.percentile(largest_losses, (1 - 1/250) * 100)

которое в настоящее время неточно и составляет 5000, поэтому нужен процесс, который может выполняться за несколько минут или меньше в течение ~ 100 000 лет.

Подробнее здесь: https://stackoverflow.com/questions/790 ... sands-of-i