用Python模拟人寿保险健康状态转移：从马尔可夫链到稳态预测（附完整代码）

用Python模拟人寿保险健康状态转移从马尔可夫链到稳态预测附完整代码在保险精算和风险管理领域预测投保人群体的健康状态演变是一项基础而关键的工作。想象你是一位保险公司的数据科学家管理层需要你评估某款寿险产品的长期风险——这本质上是一个关于人群健康状态如何随时间迁移的概率问题。而马尔可夫链Markov Chain这个以俄罗斯数学家安德雷·马尔可夫命名的概率模型恰好提供了完美的数学工具来描述这种状态转移过程。与传统统计分析不同马尔可夫链的魅力在于它的无记忆性——下一个状态只取决于当前状态与历史路径无关。这种特性使得复杂的状态转移问题变得可建模、可计算。本文将带你用Python从零构建健康状态预测模型通过代码实现两状态健康/疾病和三状态健康/疾病/死亡的马尔可夫模拟最终计算出无论初始状态如何都会收敛的稳态分布。这些技术同样适用于用户行为分析、金融信用评级等场景。1. 环境准备与基础概念在开始编写代码前我们需要明确几个核心概念。马尔可夫链由三个要素构成状态空间所有可能状态的集合如健康、疾病转移矩阵描述状态间转移概率的方阵初始分布各个状态的起始概率安装必要的Python库pip install numpy pandas matplotlib先看一个最简单的健康-疾病两状态系统。假设年度转移概率如下健康→健康0.8健康→疾病0.2疾病→健康0.7疾病→疾病0.3对应的转移矩阵P为import numpy as np P np.array([ [0.8, 0.2], # 健康行的转移概率 [0.7, 0.3] # 疾病行的转移概率 ])注意转移矩阵的行表示当前状态列表示下一状态每行概率之和必须严格等于12. 两状态马尔可夫链实现2.1 状态转移模拟让我们模拟一个投保人未来5年的健康状态变化。假设初始状态为健康概率100%def simulate_markov(transition_matrix, initial_state, steps): states [initial_state] current_state initial_state for _ in range(steps): next_state np.random.choice( len(transition_matrix), ptransition_matrix[current_state] ) states.append(next_state) current_state next_state return states # 状态编码0健康1疾病 initial 0 trajectory simulate_markov(P, initial, 5) print(状态序列:, [健康 if s0 else 疾病 for s in trajectory])单次模拟可能输出状态序列: [健康, 健康, 疾病, 健康, 疾病, 健康]2.2 可视化状态转移为了更直观理解我们绘制状态转移概率图import matplotlib.pyplot as plt from matplotlib.patches import Circle, Arrow def plot_markov_chain(P): fig, ax plt.subplots(figsize(8,4)) # 绘制状态节点 health Circle((1,1), 0.2, fclightgreen, ecblack) disease Circle((3,1), 0.2, fcsalmon, ecblack) ax.add_patch(health) ax.add_patch(disease) # 添加转移箭头 arrowprops dict(arrowstyle-, connectionstylearc3) ax.annotate(, xy(1.2,1), xytext(2.8,1), arrowpropsarrowprops) ax.annotate(, xy(2.8,1), xytext(1.2,1), arrowpropsarrowprops) # 标注概率 ax.text(2, 1.1, f{P[0,1]:.2f}, hacenter) ax.text(2, 0.9, f{P[1,0]:.2f}, hacenter) ax.text(1, 1.2, f自循环 {P[0,0]:.2f}, hacenter) ax.text(3, 1.2, f自循环 {P[1,1]:.2f}, hacenter) ax.set_xlim(0,4) ax.set_ylim(0,2) ax.axis(off) plt.show() plot_markov_chain(P)2.3 计算稳态分布马尔可夫链的一个重要特性是最终会收敛到稳态分布。计算稳态分布有三种方法方法1矩阵幂次法def steady_state_power(P, tol1e-8): n 1 while True: P_new np.linalg.matrix_power(P, n) P_next np.linalg.matrix_power(P, n1) if np.allclose(P_new, P_next, atoltol): return P_new[0] n 1 steady_power steady_state_power(P) print(幂次法稳态分布:, steady_power)方法2特征向量法def steady_state_eigen(P): eigenvalues, eigenvectors np.linalg.eig(P.T) idx np.argmin(np.abs(eigenvalues - 1)) steady eigenvectors[:,idx].real return steady / steady.sum() steady_eigen steady_state_eigen(P) print(特征向量法稳态分布:, steady_eigen)方法3解线性方程组def steady_state_linear(P): n P.shape[0] A np.vstack([(P.T - np.eye(n))[:-1], np.ones(n)]) b np.zeros(n) b[-1] 1 return np.linalg.solve(A, b) steady_linear steady_state_linear(P) print(线性方程组法稳态分布:, steady_linear)三种方法都应输出近似结果[0.777..., 0.222...]表示长期来看约77.8%的时间处于健康状态22.2%处于疾病状态。3. 三状态模型含吸收态现在扩展模型到更现实的三种状态健康(0)、疾病(1)、死亡(2)。死亡是吸收态——一旦进入就无法离开。假设转移矩阵为P_3state np.array([ [0.8, 0.18, 0.02], # 健康 [0.25, 0.65, 0.1], # 疾病 [0, 0, 1] # 死亡(吸收态) ])3.1 吸收链分析对于含吸收态的马尔可夫链我们关心两个核心问题进入吸收态的平均时间最终被吸收的概率计算基本矩阵N和吸收概率BQ P_3state[:2, :2] # 非吸收态部分 I np.eye(2) N np.linalg.inv(I - Q) # 基本矩阵 print(平均转移次数:\n, N) R P_3state[:2, 2:] B N R # 吸收概率 print(\n吸收概率:\n, B)输出示例平均转移次数: [[5.26315789 1.57894737] [2.63157895 3.15789474]] 吸收概率: [[1.] [1.]]表示从健康状态出发平均需要5.261.58≈6.84次转移被吸收从疾病状态出发平均需要2.633.16≈5.79次转移。两种状态最终都会被死亡态吸收概率100%。3.2 可视化多步转移观察不同初始条件下状态概率的演变def multi_step_probs(P, initial, steps): probs [initial] for _ in range(steps): probs.append(probs[-1] P) return np.array(probs) initial_healthy np.array([1, 0, 0]) probs multi_step_probs(P_3state, initial_healthy, 50) plt.figure(figsize(10,6)) plt.plot(probs[:,0], label健康) plt.plot(probs[:,1], label疾病) plt.plot(probs[:,2], label死亡) plt.xlabel(时间步) plt.ylabel(概率) plt.legend() plt.show()图像将清晰显示死亡概率随时间单调递增最终趋近于1。4. 实际业务应用4.1 保费定价参考假设健康状态年赔付成本500元疾病状态年赔付成本5000元死亡一次性赔付100,000元计算长期期望成本steady steady_state_linear(P_3state[:2,:2]) # 忽略死亡态 cost np.dot(steady, [500, 5000]) print(年期望赔付成本:, cost) death_prob P_3state[0,2] P_3state[1,2]*steady[1] print(年死亡率:, death_prob)4.2 敏感性分析转移概率的微小变化会如何影响稳态我们测试健康保持概率从0.7到0.9的影响health_probs np.linspace(0.7, 0.9, 5) results [] for p in health_probs: P_temp np.array([ [p, 0.3-p/10, 0.02], [0.25, 0.65, 0.1], [0, 0, 1] ]) steady steady_state_linear(P_temp[:2,:2]) results.append(steady[0]) plt.plot(health_probs, results) plt.xlabel(健康保持概率) plt.ylabel(稳态健康概率) plt.show()4.3 模型验证方法为确保模型准确性建议历史数据回溯测试比较模型预测与实际群体健康状态分布计算Kolmogorov-Smirnov统计量蒙特卡洛模拟def monte_carlo_sim(P, initial, trials, steps): counts np.zeros(P.shape[0]) for _ in range(trials): state initial for _ in range(steps): state np.random.choice(len(P), pP[state]) counts[state] 1 return counts / trials print(模拟结果:, monte_carlo_sim(P, 0, 10000, 100))完整代码已封装为Python类class HealthMarkovModel: def __init__(self, transition_matrix): self.P np.array(transition_matrix) self.n_states self.P.shape[0] def simulate(self, initial_state, steps): states [initial_state] current initial_state for _ in range(steps): current np.random.choice(self.n_states, pself.P[current]) states.append(current) return states def steady_state(self, methodeigen): if method eigen: vals, vecs np.linalg.eig(self.P.T) vec vecs[:, np.argmin(np.abs(vals - 1))].real return vec / vec.sum() elif method power: n 1 while True: Pn np.linalg.matrix_power(self.P, n) Pn1 np.linalg.matrix_power(self.P, n1) if np.allclose(Pn, Pn1, atol1e-8): return Pn[0] n 1 def absorption_analysis(self): if not np.any(np.diag(self.P) 1): raise ValueError(No absorbing states detected) absorb_mask np.diag(self.P) 1 Q self.P[~absorb_mask, :][:, ~absorb_mask] R self.P[~absorb_mask, :][:, absorb_mask] N np.linalg.inv(np.eye(Q.shape[0]) - Q) B N R return N, B