用Python+NumPy手把手模拟人寿保险健康状态预测（附完整代码）

马尔可夫链作为描述随机过程的数学工具，在保险精算领域有着广泛的应用价值。想象一下，当你需要预测未来几年内投保人群的健康状态分布时，马尔可夫链能提供清晰的数学模型。本文将带你用Python和NumPy库，从零开始构建一个人寿保险健康状态预测模型，涵盖两状态（健康/疾病）和三状态（健康/疾病/死亡）的完整实现。

1. 环境准备与基础概念

在开始编码前，我们需要确保环境配置正确。建议使用Python 3.8+版本，并安装以下依赖库：

pip install numpy matplotlib

马尔可夫链的核心特征是 无记忆性 ——下一状态的概率分布仅取决于当前状态。在人寿保险场景中，这意味着投保人明年的健康状态只与今年的状态相关，而与更早的历史无关。

关键数学概念 ：

• 状态向量：表示当前各状态概率的向量，如 [0.8, 0.2] 表示80%概率健康，20%概率患病
• 转移矩阵：描述状态间转换概率的方阵，行表示当前状态，列表示下一状态

2. 两状态模型实现

我们先实现健康-疾病两状态模型。假设转移概率如下：

• 健康→健康：0.8
• 健康→疾病：0.2
• 疾病→健康：0.7
• 疾病→疾病：0.3

对应的转移矩阵为：

import numpy as np

P = np.array([
    [0.8, 0.2],  # 健康→健康, 健康→疾病
    [0.7, 0.3]   # 疾病→健康, 疾病→疾病
])

2.1 状态预测函数

编写核心的预测函数，计算n年后的状态分布：

def predict_state(initial_state, P, n):
    """
    initial_state: 初始状态向量
    P: 转移矩阵
    n: 预测年数
    """
    current_state = initial_state.copy()
    for _ in range(n):
        current_state = np.dot(current_state, P)
    return current_state

2.2 不同初始状态的演化

我们比较三种初始状态的长期演变：

# 三种初始状态
initial_healthy = np.array([1, 0])      # 初始健康
initial_sick = np.array([0, 1])         # 初始患病
initial_mixed = np.array([0.75, 0.25])  # 混合状态

# 计算5年后的状态
years = 5
print("初始健康:", predict_state(initial_healthy, P, years))
print("初始患病:", predict_state(initial_sick, P, years))
print("混合初始:", predict_state(initial_mixed, P, years))

运行结果会显示，无论初始状态如何，长期都会收敛到 [0.777..., 0.222...] ，这与理论分析的正则马尔可夫链性质一致。

2.3 可视化演化过程

使用Matplotlib绘制状态变化曲线：

import matplotlib.pyplot as plt

def plot_evolution(initial_state, P, years, title):
    healthy_probs = []
    sick_probs = []
    current = initial_state
    for _ in range(years+1):
        healthy_probs.append(current[0])
        sick_probs.append(current[1])
        current = np.dot(current, P)
    
    plt.figure(figsize=(10,6))
    plt.plot(healthy_probs, label='健康概率', marker='o')
    plt.plot(sick_probs, label='疾病概率', marker='s')
    plt.title(title)
    plt.xlabel('年数')
    plt.ylabel('概率')
    plt.legend()
    plt.grid(True)
    plt.show()

# 绘制三种初始状态的演化
plot_evolution(initial_healthy, P, 10, "初始健康状态下的演化")
plot_evolution(initial_sick, P, 10, "初始疾病状态下的演化")
plot_evolution(initial_mixed, P, 10, "混合初始状态下的演化")

3. 三状态模型进阶

现实中人寿保险需要考虑死亡状态。扩展为三状态模型：

• 健康→健康：0.8
• 健康→疾病：0.18
• 健康→死亡：0.02
• 疾病→健康：0.25
• 疾病→疾病：0.65
• 疾病→死亡：0.1
• 死亡→死亡：1（吸收态）

转移矩阵变为：

P_3state = np.array([
    [0.8, 0.18, 0.02],  # 健康
    [0.25, 0.65, 0.1],  # 疾病
    [0, 0, 1]           # 死亡(吸收态)
])

3.1 吸收态处理技巧

死亡作为吸收态，会导致矩阵幂运算效率低下。我们可以利用分块矩阵性质优化计算：

def predict_absorbing(initial_state, P, max_years):
    """
    处理含吸收态的长期预测
    """
    # 分离非吸收态和吸收态
    Q = P[:2, :2]  # 非吸收态部分
    R = P[:2, 2:]  # 到吸收态部分
    
    # 计算基本矩阵
    N = np.linalg.inv(np.eye(2) - Q)
    
    # 计算吸收概率
    B = np.dot(N, R)
    
    results = []
    current = initial_state[:2]
    for year in range(max_years+1):
        # 计算当前非吸收态概率
        temp = np.dot(current, np.linalg.matrix_power(Q, year))
        # 计算死亡概率
        death_prob = 1 - temp.sum()
        results.append(np.append(temp, death_prob))
    
    return np.array(results)

3.2 长期趋势分析

运行三状态模型预测：

initial_healthy_3 = np.array([1, 0, 0])
initial_sick_3 = np.array([0, 1, 0])
initial_mixed_3 = np.array([0.75, 0.25, 0])

results_healthy = predict_absorbing(initial_healthy_3, P_3state, 50)
results_sick = predict_absorbing(initial_sick_3, P_3state, 50)
results_mixed = predict_absorbing(initial_mixed_3, P_3state, 50)

可视化三状态演化：

def plot_3state(results, title):
    plt.figure(figsize=(12,7))
    years = np.arange(len(results))
    plt.plot(years, results[:,0], label='健康', linestyle='-')
    plt.plot(years, results[:,1], label='疾病', linestyle='--')
    plt.plot(years, results[:,2], label='死亡', linestyle=':')
    plt.title(title)
    plt.xlabel('年数')
    plt.ylabel('概率')
    plt.legend()
    plt.grid(True)
    plt.show()

plot_3state(results_healthy, "初始健康的三状态演化")
plot_3state(results_sick, "初始疾病的三状态演化")
plot_3state(results_mixed, "混合初始的三状态演化")

4. 实用技巧与性能优化

4.1 矩阵运算加速

对于大规模状态空间，可以使用稀疏矩阵：

from scipy.sparse import csr_matrix

P_sparse = csr_matrix(P_3state)
# 稀疏矩阵乘法
result = initial_healthy_3.dot(P_sparse.power(10))

4.2 常见错误排查

维度不匹配错误 ：

# 错误示例：状态向量形状不正确
wrong_state = np.array([1, 0, 0])  # 形状(3,)
P = np.array([[0.8,0.2],[0.7,0.3]])  # 形状(2,2)
# 正确做法：确保形状匹配
correct_state = np.array([[1, 0]])  # 形状(1,2)

概率不守恒检查 ：

def validate_matrix(P):
    # 检查每行概率和是否为1
    if not np.allclose(P.sum(axis=1), 1):
        raise ValueError("转移矩阵行概率和必须为1")

4.3 实际应用扩展

将模型封装为保险产品评估类：

class InsuranceModel:
    def __init__(self, states, transition_matrix):
        self.states = states
        self.P = np.array(transition_matrix)
        validate_matrix(self.P)
        
    def simulate(self, initial_dist, years):
        """模拟多年演变"""
        if len(initial_dist) != len(self.states):
            raise ValueError("初始分布与状态数不匹配")
            
        results = np.zeros((years+1, len(self.states)))
        results[0] = initial_dist
        for year in range(1, years+1):
            results[year] = np.dot(results[year-1], self.P)
        return results

# 使用示例
states = ["健康", "疾病", "死亡"]
model = InsuranceModel(states, P_3state)
results = model.simulate([1,0,0], 30)

5. 完整代码整合

以下是可直接运行的Jupyter Notebook完整代码：

# 人寿保险马尔可夫模型完整实现
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix

class InsuranceMarkovModel:
    def __init__(self, states, transition_matrix):
        """
        初始化保险马尔可夫模型
        :param states: 状态名称列表
        :param transition_matrix: 转移概率矩阵
        """
        self.states = states
        self.P = np.array(transition_matrix)
        self._validate_matrix()
        
    def _validate_matrix(self):
        """验证转移矩阵有效性"""
        if not np.allclose(self.P.sum(axis=1), 1, atol=1e-6):
            raise ValueError("转移矩阵每行概率和必须为1")
        if self.P.shape[0] != self.P.shape[1]:
            raise ValueError("转移矩阵必须是方阵")
        if len(self.states) != self.P.shape[0]:
            raise ValueError("状态数量与转移矩阵维度不匹配")
    
    def simulate(self, initial_dist, years):
        """
        模拟状态演变
        :param initial_dist: 初始概率分布
        :param years: 模拟年数
        :return: (years+1) x n_states 的结果矩阵
        """
        if len(initial_dist) != len(self.states):
            raise ValueError("初始分布长度与状态数不匹配")
        if not np.isclose(sum(initial_dist), 1, atol=1e-6):
            raise ValueError("初始概率分布和必须为1")
            
        results = np.zeros((years+1, len(self.states)))
        results[0] = initial_dist
        for year in range(1, years+1):
            results[year] = np.dot(results[year-1], self.P)
        return results
    
    def plot_evolution(self, results, title=None):
        """绘制状态演化曲线"""
        plt.figure(figsize=(12,7))
        years = np.arange(len(results))
        for i, state in enumerate(self.states):
            plt.plot(years, results[:,i], label=state)
        
        if title:
            plt.title(title)
        plt.xlabel('年数')
        plt.ylabel('概率')
        plt.legend()
        plt.grid(True)
        plt.show()

# 两状态模型示例
print("=== 两状态健康-疾病模型 ===")
states_2 = ["健康", "疾病"]
P_2state = [[0.8, 0.2], [0.7, 0.3]]
model_2state = InsuranceMarkovModel(states_2, P_2state)

# 模拟不同初始条件
initial_healthy = [1, 0]
initial_sick = [0, 1]
initial_mixed = [0.75, 0.25]

results_healthy = model_2state.simulate(initial_healthy, 10)
results_sick = model_2state.simulate(initial_sick, 10)
results_mixed = model_2state.simulate(initial_mixed, 10)

model_2state.plot_evolution(results_healthy, "初始健康状态演化")
model_2state.plot_evolution(results_sick, "初始疾病状态演化")
model_2state.plot_evolution(results_mixed, "混合初始状态演化")

# 三状态模型示例
print("\n=== 三状态健康-疾病-死亡模型 ===")
states_3 = ["健康", "疾病", "死亡"]
P_3state = [
    [0.8, 0.18, 0.02],
    [0.25, 0.65, 0.1],
    [0, 0, 1]
]
model_3state = InsuranceMarkovModel(states_3, P_3state)

# 模拟三状态
initial_healthy_3 = [1, 0, 0]
initial_sick_3 = [0, 1, 0]
initial_mixed_3 = [0.75, 0.25, 0]

results_healthy_3 = model_3state.simulate(initial_healthy_3, 50)
results_sick_3 = model_3state.simulate(initial_sick_3, 50)
results_mixed_3 = model_3state.simulate(initial_mixed_3, 50)

model_3state.plot_evolution(results_healthy_3, "初始健康的三状态演化")
model_3state.plot_evolution(results_sick_3, "初始疾病的三状态演化")
model_3state.plot_evolution(results_mixed_3, "混合初始的三状态演化")