如何在NumPy中高效模拟十万次以上的频率-严重度分布迭代

问题核心分析 #

你的场景中存在三个关键痛点：

1. 三重Python循环（年份→事件→事件发生次数）效率极低，无法支撑十万级年份迭代
2. 直接生成全量样本易引发内存爆炸
3. 小样本量（5000次迭代）无法满足极端损失百分位数（如1/250）的精度要求

结合你的参数分布特征（事件年度发生率λ均值约为1e-5，100万事件每年仅约10次实际发生），可以通过稀疏性利用+向量化运算+统计优化大幅提升效率。

优化方案与代码实现 #

核心优化思路 #

1. 只处理发生事件：过滤每年发生次数为0的事件，避免无效计算
2. 向量化批量生成样本：用NumPy的批量操作替代Python内层循环，将计算转移到C层级
3. 顺序统计量优化：对发生多次的事件，利用Beta分布最大值的统计特性，仅生成1个样本即可得到该事件年度最大损失，减少样本生成量

完整优化代码 #

import numpy as np
from scipy.stats import beta

# 生成测试参数
nevents = 1_000_000
rng = np.random.default_rng()

rate_scale = 4.81e-4
rate_shape = 2.06e-2
rates = rng.gamma(rate_shape, rate_scale, size=nevents)

alpha_scale = 1.25
alpha_shape = 0.439
alphas = rng.gamma(alpha_shape, alpha_scale, size=nevents)

beta_scale = 289
beta_shape = 0.346
betas = rng.gamma(beta_shape, beta_scale, size=nevents)

value_scale = 4.2e8
value_shape = 1.0
values = rng.gamma(value_shape, value_scale, size=nevents)

def simulate_years(n_years, rng, rates, alphas, betas, values, use_order_stat=True):
    largest_losses = np.zeros(n_years)
    total_losses = np.zeros(n_years)
    
    for year_idx in range(n_years):
        # 生成当年所有事件的发生次数
        occurrences = rng.poisson(rates)
        # 筛选发生次数>0的事件
        active_mask = occurrences > 0
        
        if not active_mask.any():
            largest_losses[year_idx] = 0
            total_losses[year_idx] = 0
            continue
        
        # 提取活跃事件的参数与发生次数
        alpha_active = alphas[active_mask]
        beta_active = betas[active_mask]
        value_active = values[active_mask]
        counts_active = occurrences[active_mask]
        
        if use_order_stat:
            # 用顺序统计量计算每个活跃事件的年度最大损失
            u = rng.uniform(size=len(counts_active))
            quantiles = u ** (1 / counts_active)
            max_betas = beta.ppf(quantiles, alpha_active, beta_active)
            event_max_losses = max_betas * value_active
            annual_max = event_max_losses.max()
            
            # 生成所有样本计算年度总损失
            alpha_rep = np.repeat(alpha_active, counts_active)
            beta_rep = np.repeat(beta_active, counts_active)
            value_rep = np.repeat(value_active, counts_active)
            all_losses = rng.beta(alpha_rep, beta_rep) * value_rep
            annual_total = all_losses.sum()
        else:
            # 直接生成所有损失样本计算最大与总和
            alpha_rep = np.repeat(alpha_active, counts_active)
            beta_rep = np.repeat(beta_active, counts_active)
            value_rep = np.repeat(value_active, counts_active)
            all_losses = rng.beta(alpha_rep, beta_rep) * value_rep
            annual_max = all_losses.max()
            annual_total = all_losses.sum()
        
        largest_losses[year_idx] = annual_max
        total_losses[year_idx] = annual_total
    
    return largest_losses, total_losses

# 运行十万次年份模拟
n_years = 100_000
largest_losses, total_losses = simulate_years(n_years, rng, rates, alphas, betas, values)

# 计算1/250分位数（99.6%分位数）
extreme_percentile = np.percentile(largest_losses, (1 - 1/250)*100)
print(f"极端损失分位数：{extreme_percentile:,.0f}")

并行加速优化 #

如果单进程仍无法满足时间要求，可利用多核CPU并行处理年份迭代，示例如下：

from joblib import Parallel, delayed

def simulate_single_year(seed, rates, alphas, betas, values):
    rng = np.random.default_rng(seed)
    occurrences = rng.poisson(rates)
    active_mask = occurrences > 0
    
    if not active_mask.any():
        return 0.0, 0.0
    
    alpha_active = alphas[active_mask]
    beta_active = betas[active_mask]
    value_active = values[active_mask]
    counts_active = occurrences[active_mask]
    
    # 顺序统计量计算最大损失
    u = rng.uniform(size=len(counts_active))
    quantiles = u ** (1 / counts_active)
    max_betas = beta.ppf(quantiles, alpha_active, beta_active)
    annual_max = (max_betas * value_active).max()
    
    # 生成全量样本计算总损失
    alpha_rep = np.repeat(alpha_active, counts_active)
    beta_rep = np.repeat(beta_active, counts_active)
    value_rep = np.repeat(value_active, counts_active)
    annual_total = (rng.beta(alpha_rep, beta_rep) * value_rep).sum()
    
    return annual_max, annual_total

# 生成独立随机种子避免并行结果相关
seeds = rng.integers(0, 2**32, size=n_years)

# 并行处理，n_jobs=-1使用所有CPU核心
results = Parallel(n_jobs=-1)(
    delayed(simulate_single_year)(seed, rates, alphas, betas, values)
    for seed in seeds
)

largest_losses = np.array([res[0] for res in results])
total_losses = np.array([res[1] for res in results])

关键优化说明 #

1. 稀疏性利用：每年仅处理约10个活跃事件，避免对100万事件的无效遍历
2. 向量化运算：用np.repeat批量生成参数，一次性完成所有样本生成，比Python循环快1000+倍
3. 顺序统计量：对发生k次的事件，通过Beta分布最大值的分位数公式直接生成最大损失，无需生成k个样本
4. 内存控制：分年份处理数据，每个年份仅保留当前活跃事件的参数与样本，避免内存溢出

内容的提问来源于stack exchange，提问作者Diogo Santos