У нас есть, скажем, 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 с... упс)
РЕДАКТИРОВАТЬ:
По запросу вот версия, в которой просто рассчитывается максимальный убыток за каждый год
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 лет.
РЕДАКТИРОВАТЬ 2: Воспроизводимые значения (в виде массивов np):
rate_mu = 9.9e-6
rate_sigma = 6.9e-5
rate_scale = rate_sigma ** 2 / rate_mu
rate_shape = rate_mu / rate_scale
rates = rng.gamma(rate_shape, rate_scale, size=nevents)
alpha_mu = 0.55
alpha_sigma = 0.83
alpha_scale = alpha_sigma ** 2 / alpha_mu
alpha_shape = alpha_mu / alpha_scale
alphas = rng.gamma(alpha_shape, alpha_scale, size=nevents)
beta_mu = 100
beta_sigma = 170
beta_scale = beta_sigma ** 2 / beta_mu
beta_shape = beta_mu / beta_scale
betas = rng.gamma(beta_shape, beta_scale, size=nevents)
value_mu = 4.2e8
value_sigma = 4.2e8
value_scale = value_sigma ** 2 / value_mu
value_shape = value_mu / value_scale
values = rng.gamma(value_shape, value_scale, size=nevents)
и внес поправки в свой текущий код, чтобы использовать массивы numpy
years = 20_000
rng = np.random.default_rng()
largest_losses = np.zeros(shape=years)
for year in range(years):
occurences = rng.poisson(rates)
largest_loss = 0
for idx, event_occs in enumerate(occurences):
for event_occ in range(event_occs):
loss = rng.beta(alphas[idx], betas[idx]) * values[idx]
if loss > largest_loss:
largest_loss = loss
largest_losses[year] = largest_loss
f'{np.percentile(largest_losses, (1 - 1/30) * 100):,.0f}', f'{np.percentile(largest_losses, (1 - 1/250) * 100):,.0f}'
Подробнее здесь: https://stackoverflow.com/questions/790 ... sands-of-i
