简介：2024年美国大学生数学建模竞赛ICM Problem E赛题的Python实现方案，聚焦财产保险可持续性建模，面向数学建模参赛者、Python数据分析入门或进阶学习者，也适用于毕业设计、课程设计、大作业和工程实训场景。资源包共43个文件，大小18.43MB，包含7个Python脚本（如灰色预测、模糊综合评价、层次分析、支持向量机等模型实现）、3个Excel数据集、2个PDF赛题文档、多张可视化PNG/SVG图表及PSD设计源文件，代码、数据与文档组织清晰。已有111人学习下载。方案完整覆盖从数据读取、模型训练到灵敏度分析、结果可视化的流程，可直接运行学习；图表可复用于论文或答辩展示，适合作为建模竞赛复盘或项目起步的参考。

1. 财产保险可持续性：美赛这道题在考什么，Python建模从哪下手

每年美赛的财产保险可持续性题都会劝退一批习惯跑机器学习的人——他们翻遍数据找标签和准确率，最后发现题目问的根本不是“哪张保单会出险”，而是“一家财产保险公司在现有保费、赔付与资本结构下，还能不能继续经营下去”。用 Python 做这道题，核心不是堆模型，而是把产险精算里那套现金流逻辑搬进蒙特卡洛仿真：从保单的理赔频率和强度分布出发，模拟未来多年的赔付总额，计算破产概率与偿付能力充足率。适合两类读者：正在备赛、想要可复现 baseline 的数模学生，以及想用 Python 验证保险业务可持续性的从业者。这篇按“拆问题→拟合分布→现金流仿真→避坑→敏感性分析”的顺序展开，中间每一步都给可跑的代码和参数说明。

2. 先拆可持续性：保费、赔付与偿付能力怎么定量

“可持续”这个词在数模题里最容易被写成口号。要让它能被 Python 计算，第一步是把概念翻译成三个可量化的口径，论文里所有结论都挂在它们下面。

2.1 三个量化口径：承保亏损、偿付能力与破产概率

第一个口径是承保结果。财产保险公司每年收上来一批保费，赔出去一批赔款，二者之比叫损失率（loss ratio），加上费用率之后再叫综合成本率。损失率超过 100%，说明承保端在亏钱，只能靠投资收益弥补。这个口径在 Python 里就是一行除法，但它决定了后续所有模拟的方向：保费不是拍脑袋定的，而是要和赔付分布对齐。

第二个口径是偿付能力充足率。监管会要求保险公司保留一笔最低资本，用于吸收超出预期的赔付。用模拟盈余与要求资本做比较，盈余低于要求的次数占比，就是偿付能力不足的概率。这个口径比损失率更接近“可持续”，因为它把时间跨度和波动性都包含了进来。

第三个口径是长期生存概率，也叫破产概率（ruin probability）。给定初始资本、每年保费、理赔频率与强度分布，模拟未来 10 到 20 年现金流，记录盈余第一次小于零的时间。这是蒙特卡洛最擅长回答的问题，也是美赛评阅人最容易看懂的量化结果。

口径 | 定义 | Python 里怎么算 | 适合回答的问题 承保损失率 | 总赔付 / 总保费 | 按年聚合后直接除 | 现在的费率够不够覆盖理赔 偿付能力充足率 | 模拟盈余 / 最低要求资本 | 蒙特卡洛每年计算比率 | 资本缓冲是否充足 破产概率 | 盈余首次为负的模拟占比 | 多年现金流循环统计 | 长期经营是否可持续

三个口径不是互相替代，而是从近到远的三个观察距离。竞赛时间有限，我会建议把重心放在破产概率上，损失率和偿付能力比率作为辅助指标并列给出。

2.2 为什么总赔付要拆成“频率×强度”建模

很多第一次做题的人会把各年总赔付直接当成一个分布去拟合，然后用它模拟未来。这种做法在数据量大的时候看起来可行，但有两个问题：第一，总赔付把不同风险特征的保单混在一起，方差被严重低估；第二，它没法回答“如果费率上调 20%，破产概率降多少”这类问题，因为费率调整改变的是承保门槛和保单结构，而不是总赔付的统计形状。

更稳妥的做法是把每张保单的年赔付拆成两部分：赔付次数（频率）和单次赔付额（强度）。每张保单一年里可能赔 0 次、1 次或多次，赔一次的钱也可能从几千到几十万不等。总赔付就写成复合分布的形式：

S = Σᵢ Σⱼ₌₁ᴺⁱ Xᵢⱼ

其中 Nᵢ 是第 i 张保单在一年里的理赔次数，Xᵢⱼ 是第 i 张保单第 j 次理赔的金额。Nᵢ 通常是泊松分布或负二项分布，Xᵢⱼ 常用对数正态、Gamma 或 Weibull。把频率和强度分开之后，既能体现不同地区、不同房屋类型的风险差异，也能在模拟时为每张保单独立抽样，尾部风险比直接拟合计总额更真实。

2.3 选型理由：泊松-对数正态作为 baseline

我一般会把“泊松频率 + 对数正态强度”作为 baseline 组合，理由有三个。其一，泊松分布只有一个参数 λ，可以直接用保单特征做回归扩展成 λᵢ，解释性很强；其二，对数正态分布在保险精算里用得最广，能捕捉“多数赔付很少、少数赔付很大”的右偏形态，参数也只有两个；其三，scipy 和 statsmodels 对这两个分布的支持最成熟，拟合和抽样都是一行代码的事。

备选方案是负二项分布替代泊松。负二项多一个离散参数，能处理“实际数据比泊松预测更分散”的情况（也叫过度离散），赛题数据如果理赔次数方差远大于均值，就用负二项。强度方面，Gamma 分布尾部比对数正态更薄，适合赔付上限明显、大额极端值较少的场景。选哪个不是越复杂越好，而是要看拟合优度检验结果，下一章直接演示判断方法。

3. Python 实现损失分布拟合：从保单数据到赔付分布参数

题目给的原始数据通常不是现成的统计样本，而是理赔明细和保单明细两份表。想把分布参数估准，得先把数据整理成“每张保单在观测期内的理赔次数和理赔总额”。

3.1 数据预处理：把理赔记录归位到保单年度

常见的赛题数据结构是：保单表含保单编号、起保日期、保险金额、房屋类型、地区等字段；理赔表含保单编号、理赔日期、赔付金额。要做的是把每张理赔单匹配到它所属的保单年度里，再按保单聚合出理赔次数和赔付总额。

import pandas as pd import numpy as np # 保单表和理赔表读入后统一日期格式 policy_df["policy_start"] = pd.to_datetime(policy_df["policy_start"]) policy_df["policy_end"] = pd.to_datetime(policy_df["policy_end"]) claim_df["claim_date"] = pd.to_datetime(claim_df["claim_date"]) # 判断理赔是否发生在保单有效期内 merged = claim_df.merge(policy_df, on="policy_id", how="left") merged["is_covered"] = ( (merged["claim_date"] >= merged["policy_start"]) & (merged["claim_date"] <= merged["policy_end"]) ) # 只保留有效理赔，按保单汇总 valid_claims = merged[merged["is_covered"]] policy_stats = valid_claims.groupby("policy_id").agg( claim_count=("claim_amount", "count"), claim_sum=("claim_amount", "sum") ).reset_index() # 和保单主表拼接，保单主表决定每一“张”样本 df_policy = policy_df.merge(policy_stats, on="policy_id", how="left") df_policy["claim_count"] = df_policy["claim_count"].fillna(0) df_policy["claim_sum"] = df_policy["claim_sum"].fillna(0)

这段代码做了三件事：日期标准化、理赔归属判断、按保单聚合。关键是is_covered这一步，如果理赔日期恰好落在保单起止时间之外还在统计口径里，后面拟合出来的频率参数就会虚高。另外注意how="left"保留了没有理赔的保单，这些样本对频率估计非常重要——它们代表“一年没出险”的绝大多数情况，不能因为没赔钱就被丢进黑匣子。

3.2 频率与强度分别拟合：泊松回归与对数正态

理赔次数可以直接用全部保单样本拟合。赛题数据里保单带有地区、房屋类型、保险金额等特征，用泊松回归比用单一 λ 更合理。理赔强度则只用“理赔次数大于 0”的样本拟合，因为零赔付样本对单次赔付金额的分布没有信息。

import statsmodels.api as sm import statsmodels.formula.api as smf from scipy import stats # 频率模型：理赔次数 ~ 特征，选择 Poisson 连接 formula_freq = "claim_count ~ region + building_type + np.log(insured_value)" model_freq = smf.glm( formula=formula_freq, data=df_policy, family=sm.families.Poisson() ).fit() # 提取强度样本：只保留有赔付记录的保单 severity_data = valid_claims["claim_amount"].values # 拟合对数正态分布：floc=0 强制位置参数为 0 shape, loc, scale = stats.lognorm.fit(severity_data, floc=0) sigma = shape # 对数标准差 mu = np.log(scale) # 对数均值 print(f"泊松回归 AIC: {model_freq.aic:.2f}") print(f"对数正态参数: mu={mu:.4f}, sigma={sigma:.4f}")

参数说明：smf.glm的formula是 statsmodels 的类 R 语法，np.log(insured_value)表示对保险金额取对数后作为线性项，这符合“保额越高理赔可能越多”的实际关系；family=sm.families.Poisson()指定泊松分布，默认连接函数是 log。强度部分，stats.lognorm.fit返回三个参数，shape是对数标准差 σ，scale等于 e^μ，所以取对数就得到 μ。这里最常见的坑是忘了floc=0，如果放任 loc 自由估计，小赔付样本会被一个偏移量强行吸收，后续模拟出来的赔付金额下限会变成负数。

如果理赔次数的方差明显大于均值，把family换成sm.families.NegativeBinomial()再拟一次，比较 AIC 谁小用谁。AIC 是赤池信息准则，用于衡量拟合与复杂度之间的平衡，不是看单个参数大小，而是看不同模型之间的相对值。

3.3 拟合优度检验：KS 检验与 Q-Q 图

分布选型不能靠肉眼，我用两个手段交叉验证。KS 检验给出一个 p 值，p 大于 0.05 表示不能拒绝“样本来自该分布”的原假设；Q-Q 图则用来看偏离发生在哪个部位——尾部还是中间。

from scipy import stats # KS 检验：注意参数顺序与 lognorm 定义一致 ks_stat, ks_p = stats.kstest( severity_data, "lognorm", args=(sigma, 0, np.exp(mu)) ) # Q-Q 图：样本分位数 vs 理论分位数 import matplotlib.pyplot as plt theoretical = stats.lognorm.ppf( np.linspace(0.01, 0.99, len(severity_data)), sigma, 0, np.exp(mu) ) sample_sorted = np.sort(severity_data) plt.figure(figsize=(6, 5)) plt.scatter(theoretical, sample_sorted, s=8) plt.plot([sample_sorted.min(), sample_sorted.max()], [sample_sorted.min(), sample_sorted.max()], "r--", lw=1) plt.xlabel("理论分位数") plt.ylabel("样本分位数") plt.title("Q-Q Plot: 对数正态拟合") plt.show()

kstest的args必须按照(shape, loc, scale)的顺序传，和lognorm.fit的返回顺序一致，写反了 p 值会变得毫无意义。Q-Q 图在高分位处如果明显偏离直线，说明尾部被低估，这时候要换 Gamma 或 Burr 分布重新拟合。拟合优度检验不是一次性通过就结束，它是整个建模流程的守门员，后面所有仿真结论都建立在这些分布“足够像真实数据”的前提上。

4. 用蒙特卡洛仿真估算持续经营的概率：现金流模拟代码

分布参数拟合好之后，进入整篇论文最核心也最出数的地方：模拟未来多年的现金流，统计破产概率。这个环节的代码量不大，但每个参数都直接影响结论，需要解释清楚。

4.1 一年总赔付的复合分布抽样

每年总赔付的模拟分两层：先给每张保单抽样出理赔次数，再给每一次理赔抽样出赔付金额。这样保单之间的风险差异自然就保留了下来，不需要人为指定一个总赔付方差。

# 从拟合好的模型中生成一组“未来保单” n_policy = df_policy.shape[0] def simulate_one_year(policy_data, mu, sigma, seed=0): rng = np.random.default_rng(seed) # 频率：根据保单特征预测泊松强度 lam = model_freq.predict(policy_data) freq = rng.poisson(lam) # 强度：对每个理赔抽一次对数正态 total_loss = 0.0 n_claims = freq.sum() if n_claims > 0: severities = rng.lognormal(mean=mu, sigma=sigma, size=n_claims) total_loss = severities.sum() return total_loss, n_claims

这里有一个很容易忽略的地方：泊松回归的predict返回的是每张保单的平均理赔次数 λ，不是次数本身。用rng.poisson(lam)才能生成“这张保单今年到底赔几次”的随机数。理赔强度抽样用了rng.lognormal，注意它的mean和sigma是正态分布参数，不是对数正态的均值，所以直接传第 3 章拟合出来的mu和sigma就是对的。如果赛题里有免赔额和保额上限，强度抽样之后还要做截断处理：赔付额 = min(max(x - 免赔额, 0), 保额上限 - 免赔额)，这一步漏掉会让模拟赔付明显偏高。

4.2 多年现金流与变破产概率计算

破产概率的算法是：给一个初始盈余，假设每年保费收入固定（或用预测模型逐年调整），赔付支出按上面的模拟函数生成，加上投资收益和固定费用，更新盈余。只要某一年盈余小于零，就记一次破产，然后进入下一次模拟。

def simulate_ruin_probability(policy_data, mu, sigma, annual_premium, initial_capital, years=10, sims=2000, invest_rate=0.03, fixed_cost_rate=0.15, seed=42): rng = np.random.default_rng(seed) ruin_count = 0 ruin_year_records = [] for sim_id in range(sims): surplus = initial_capital ruined = False for year in range(years): # 当年赔付（复合分布抽样） total_loss, _ = simulate_one_year( policy_data, mu, sigma, seed=seed + sim_id * 100 + year ) # 现金流更新：先投资增值，再加保费，减赔付和费用 surplus = surplus * (1 + invest_rate) surplus += annual_premium * (1 - fixed_cost_rate) surplus -= total_loss if surplus < 0: ruined = True ruin_year_records.append(year + 1) break if ruined: ruin_count += 1 ruin_prob = ruin_count / sims return ruin_prob, ruin_year_records ruin_prob, ruin_years = simulate_ruin_probability( policy_data=df_policy, mu=mu, sigma=sigma, annual_premium=df_policy["premium"].sum(), initial_capital=df_policy["premium"].sum() * 1.2, years=10, sims=2000 ) print(f"10 年破产概率: {ruin_prob:.3f}")

参数说明：initial_capital是初始资本金，我一般取当年保费的 1.2 到 2 倍，这个比例对应监管最低资本要求；invest_rate是投资收益率，取 0.03 比较保守，想分析利率影响就把它改成可调参数；fixed_cost_rate是运营费用率，0.15 是产险业的常见区间。seed拆成seed + sim_id * 100 + year是为了保证每次模拟的保单抽样序列不重复，同时保证整体可复现。这里的annual_premium在第一个版本里设成恒定值，实际论文里可以加一个“费率逐年上涨”的分支，看看多高的涨幅能把破产概率压到可接受水平。

4.3 结果解读：把概率翻译成可持续性结论

模拟完不能只报一个破产概率了事，评阅人要看这张概率表怎么指导决策。常见补法是输出赔付分布的分位数表，把它和保费、资本并列对照。

# 只看首年赔付分布，用来评估保费的充足程度 first_year_losses = [] for i in range(5000): loss, _ = simulate_one_year(df_policy, mu, sigma, seed=10000 + i) first_year_losses.append(loss) loss_arr = np.array(first_year_losses) print(f"首年赔付均值: {loss_arr.mean():,.0f}") print(f"95%分位数: {np.quantile(loss_arr, 0.95):,.0f}") print(f"99%分位数: {np.quantile(loss_arr, 0.99):,.0f}")

如果 95% 分位数超过当年保费，就说明“保费储备不足”是有量化证据的：一百年里大概有五年会赔穿。反过来，如果 99% 分位数远低于保费，说明定价偏高，可以下调费率换续保率。这三个数字放进论文里，可持续性的结论就有了锚点。到这一步，模型主流程已经完整，但真正决定论文分数的往往是下一章那些细节坑。

5. 建模避坑：影响最终结论的四个常见问题

5.1 忽略巨灾风险的相关性，尾部被严重低估

现象：模拟出来的破产概率非常低，例如只有 0.5%，但实际数据里某一年某个地区同时发生多起大额理赔，单年赔付直接超过保费的 1.8 倍。

原因：默认了所有保单的理赔相互独立。同一地区、同一片洪水或台风影响范围内的保单，风险高度相关，用独立假设抽样等于把尾部风险人为摊薄了。

解决：在模拟中给保单引入一个区域共因变量。比如每年生成一个区域风险因子 R，同区域所有保单的 λ 同时乘上这个因子，因子本身服从对数正态分布。实现起来只需在simulate_one_year里加一步：lam = model_freq.predict(policy_data) * region_factor。这样相关结构就有了，尾部也能体现“区域性巨灾同时爆发的场景”。

5.2 把不同风险的保单混在一起做总赔付拟合

现象：KS 检验 p 值远大于 0.05，看似拟合得很好，但模拟结果却和真实年度赔付对不上，尤其是高频低损和低频高损两类保单混合时。

原因：总赔付分布是多类保单的混合，混合分布虽然整体形状能被单个分布勉强逼近，但它在模拟时无法区分“大量小额赔付”和“少量大额赔付”的比例变化。

解决：使用第 3 章的分层做法——先按地区和房屋类型分组拟合强度参数，再在仿真里按保单特征分别抽样。如果赛题数据没有给保单维度特征，那就至少按“住宅类”和“商业类”拆分，拆完再拟合，模拟结果会稳定得多。

5.3 免赔额和保额上限是模拟出来的，不是数据长这样

现象：理赔强度样本用原始赔付金额拟合，模拟时也不做截断，结果模拟的平均赔付比真实赔付高 30% 以上。

原因：保单条款通常有免赔额（比如每次理赔前 5000 元自付）和保额上限。观测到的理赔金额是截断后的结果，直接用观测值拟合，再在模拟中还原原始损失，逻辑就错了。

解决：如果数据能区分“实际损失”和“赔付金额”，优先用实际损失拟合分布，在模拟时加截断逻辑：payout = min(max(loss - deductible, 0), limit - deductible)。如果数据只有赔付金额，那就把分布拟合并保持模拟时不做二次截断，只做一次免赔调整，并在论文里明确说明数据口径。这道题的数据建模口径，审稿人看得很细。

5.4 时间口径混用，暴露期长短不齐

现象：某年破产概率比别的组高出一倍，查代码没问题，查数据发现有一批保单的起保日期是年中，观测窗口只有半年，却被按一整年计算了理赔次数。

原因：数模赛题给的保单常常是陆续生效的，不是统一的 1 月 1 日到 12 月 31 日。用自然年聚合理赔次数，有些保单暴露了大半年，有些只有两个月，频率模型的 λ 被系统性拉低或拉高。

解决：在 3.1 的预处理阶段增加一个暴露时长字段exposure = (policy_end - policy_start).days / 365.25，泊松回归中使用exposure作为偏移量：formula = "claim_count ~ region + building_type + np.log(insured_value) + np.log(exposure)"，或者直接用offset=np.log(exposure)。这样不同起保日期的保单在频率上就具备了可比性，不会因为日期切得不好得出离谱结论。

6. 进阶：费率与免赔额的敏感性分析，找到性价比最高的调整

模型能跑通之后，评阅人最想看到的是“你的方案能反哺决策”。我用得最多的进阶做法是费率与免赔额的双参数敏感性分析：固定其他参数，对保费乘数和免赔额做二维网格搜索，重算破产概率，找出“既能显著降低破产概率，又不至于让费率涨幅吓跑客户”的折中带。

# 网格搜索：保费乘数 1.0~1.7，免赔额 0 / 3000 / 8000 / 15000 premium_factors = np.arange(1.0, 1.8, 0.1) deductibles = [0, 3000, 8000, 15000] results = [] for pf in premium_factors: for ded in deductibles: # 重新拟合强度时要考虑免赔额的影响，这里简化处理 adjusted_mu = mu - np.log1p(ded / np.exp(mu) * 0.3) prob, _ = simulate_ruin_probability( policy_data=df_policy, mu=adjusted_mu, sigma=sigma, annual_premium=df_policy["premium"].sum() * pf, initial_capital=df_policy["premium"].sum() * 1.2, years=10, sims=1000, seed=2024 ) results.append((pf, ded, prob)) for pf, ded, prob in sorted(results, key=lambda x: x[2]): print(f"费率乘数 {pf:.1f}, 免赔额 {ded:>6,}: 破产概率 {prob:.3f}")

注意这段里adjusted_mu是我习惯用的简化处理：免赔额提高后，可赔的平均金额下降，用一个小系数按比例调整 μ。如果时间充裕，更严谨的做法是把截断逻辑写进simulate_one_year，每种免赔额重新走一遍赔付抽样。网格搜索做完，把结果画成二维热力图或等高线图，破产概率低于 2% 的“绿色区域”就是政策建议区间。热点词“数据规范化”这里也能用上：把保费、赔付额统一按保额做归一化后再画热力图，不同方案的对比才不会受保额规模干扰。

这类赛题做到最后，我最大的教训是随机种子和参数表一定要尽早固定在论文附录里。同一个模型，种子不一样，破产概率从 3% 跳到 6% 都有可能，评阅人复现不出来，整个模型的公信力就崩了。我现在的习惯是每次跑完网格搜索都把种子、输出、参数清单存一份 CSV，论文里贴的图和数字严格取自身份一致的某一次运行。希望这些约定能帮你在赛场上少踩几个坑，把精力留给真正有区分度的结论。

本文还有配套的精品资源，点击获取