在 numpy 中有效模拟数千次迭代中的许多频率-严重性分布

Question

我在工作中遇到如下问题：

我们有 100 万个可能的事件来定义频率-严重性分布。 对于每个事件，我们都有一个定义泊松分布的年率，以及定义贝塔分布的 alpha 和 beta 参数。目标是按照>100,000“年”的顺序进行模拟，每年被定义为获取每个事件的频率 N，并获取相对 beta 分布的 N 个样本。

对我来说，麻烦的事实是如何有效地从 Beta 分布 Beta_i 中获取 N_i ~ Poisson(lambda_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 (72s -> 66s)
• 一次在循环外计算所有年份的频率（66s -> 73s...哎呀）

理想情况下，我希望它运行得尽可能快，或者进行尽可能多的迭代，而且我之前使用 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 个事件的事件池。

为了指定目标，我想看看是否可以计算损失的百分位数，例如250 分之一的损失

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

目前的精度为 5,000，因此需要一个可以在几分钟或更短时间内运行约 100,000 年的过程。

Answer 1

这是模拟年度最大损失和总损失的一种方法。 请注意，我使用普通的 NumPy 数组而不是 Pandas。

def simulate(nyears, rate, alpha, beta, value, rng):
    max_loss = np.empty(nyears, dtype=np.float64)
    total_loss = np.zeros(nyears, dtype=np.float64)
    for year in range(nyears):
        n = rng.poisson(rate)
        k = np.nonzero(n)[0]
        a = np.repeat(alpha[k], n[k])
        b = np.repeat(beta[k], n[k])
        v = np.repeat(value[k], n[k])
        losses = rng.beta(a, b)*v
        max_loss[year] = losses.max()
        total_loss[year] = losses.sum()
    return max_loss, total_loss

它可能比你的版本快一点，但我不认为它会快一个数量级。

为了测试它，我使用了这个脚本：

# Generate some test data.
rng = np.random.default_rng(121263137472525314065)
nevents = 100604
# Other than the mean of the rate being approximately 1, there is no
# special significance to the random distributions used here.
rate = rng.gamma(2, scale=0.5, size=nevents)
alpha = rng.uniform(0.5, 1.5, size=nevents)
beta = rng.uniform(1.25, 2.0, size=nevents)
value = rng.exponential(size=nevents)

nyears = 50000
max_loss, total_loss = simulate(nyears, rate, alpha, beta, value, rng)
print(np.percentile(max_loss, (1 - 1/250) * 100))