日期
Jan 15, 2025
SEIR简单个体随机模型
相较于SIR常微分方程传播模型，SEIR模型在保留最少限度的推广的同时，考虑了个体的自然史状态和于当前状态停留的时长两个属性。该模型将传播速率中的传播速率系数拆分为接触频率与单次接触感染概率的乘积，并用于个体间接触发生感染的判定。
下载讲解视频1 模型假设与符号定义
1.1 模型假设
- 个体属性：每个个体具有自然史状态（易感、潜伏、感染、康复）和在当前状态的持续时间。
- 传播机制：传播速率系数被拆分为接触频率（每天接触人数）和单次接触感染概率。
- 随机接触：个体的每日接触对象是随机选择的。
- 转化率：潜伏期和感染期的转化率分别为σ和γ。
1.2 符号定义
- m：模拟总时长（天）
- n：总人口数量（人）
- Ts：时间步长（天）
- c：日接触率（人/天）
- q：传染概率（无量纲）
- σ：潜伏期转化率（1/天）
- γ：感染期恢复率（1/天）
1.3 数据储存
(1) 个体状态（stage）
- 定义: 每个个体的健康状态，使用整数数组存储。
- 状态编码:
- 0: 易感者（S）
- 1: 潜伏者（E）
- 2: 感染者（I）
- 3: 康复者（R）
(2) 状态持续时间（durationSinceEnteredTheStage）
- 定义: 表示个体进入当前状态后的持续时间（以天数为单位）。
- 数据类型: 数值数组，长度为总人数 ( n )。
- 用途:
- 对潜伏者（E）: 判断是否转变为感染者（达到潜伏期时长）。
- 对感染者（I）: 判断是否康复（达到感染期时长）。
(3) 汇总数据
- 定义: 记录每个时间步的群体整体状态。
- 数据类型: 数值数组，长度为模拟的时间步长 ( m )。
- 记录内容:
- S_record: 易感者数量记录。
- E_record: 潜伏者数量记录。
- I_record: 感染者数量记录。
- R_record: 康复者数量记录。
- Rt_record: 实时再生数记录。
2 个体属性
2.1 属性定义
每个个体在模型中具有以下属性：
-
自然史状态（stage）：
- 0：易感者（Susceptible, S）
- 1：潜伏者（Exposed, E）
- 2：感染者（Infectious, I）
- 3：康复者（Recovered, R）
-
当前状态的持续时间（duration）：记录个体在当前状态停留的时间（天）。
2.2 人群状态初始化
模型初始化时，所有个体默认处于易感状态，部分个体被设定为初始潜伏者。
def initialize_states(self, initial_exposed=10):
"""
初始化模型状态
Args:
initial_exposed (int): 初始潜伏者数量（人）
"""
# 初始化个体状态数组（0:S, 1:E, 2:I, 3:R）
self.stage = np.zeros(self.n, dtype=int)
self.stage[0:initial_exposed] = 1 # 设置初始潜伏者
self.duration = np.zeros(self.n) # 记录个体在当前状态的持续时间（天）
# 初始化时间序列记录数组
self.S_record = np.zeros(self.m) # 易感者数量记录（人）
self.E_record = np.zeros(self.m) # 潜伏者数量记录（人）
self.I_record = np.zeros(self.m) # 感染者数量记录（人）
self.R_record = np.zeros(self.m) # 康复者数量记录（人）
self.Rt_record = np.zeros(self.m) # 实时再生数记录（无量纲）
# 新增记录数组
self.new_exposed_record = np.zeros(self.m) # 新增潜伏者数量记录（人/天）
self.new_recovered_record = np.zeros(self.m) # 新增康复者数量记录（人/天）
3 个体行为
3.1 接触传播
在每个时间步，易感者与感染者的接触会导致潜伏者的产生。具体步骤如下：
- 接触次数生成：每个易感者当天的接触次数服从泊松分布，参数为
c * ts。 - 接触者选择：随机选择接触的个体。
- 感染判断：若接触到感染者，计算感染概率
inf_prob = 1 - (1 - q) ** n_infected_contacts，并根据该概率决定是否转变为潜伏者。
def simulate(self):
"""运行SEIR模型仿真"""
for i in range(self.m):
# 记录本时间步开始时的状态，用于计算新增人数
prev_recovered = np.sum(self.stage == 3)
# 获取当前易感者和感染者的索引
susceptible_indices = np.where(self.stage == 0)[0]
infected_indices = np.where(self.stage == 2)[0]
# 记录新增潜伏者的计数器
new_exposed_count = 0
if len(infected_indices) > 0: # 只在有感染者时进行传播计算
# 为每个易感者生成接触者
for s_idx in susceptible_indices:
# 生成该易感者当天接触的人数（泊松分布）
n_contacts = np.random.poisson(self.c * self.ts)
if n_contacts > 0:
# 随机选择接触者
contacts = np.random.choice(self.n, size=n_contacts)
# 计算接触的感染者数量
n_infected_contacts = np.sum(np.isin(contacts, infected_indices))
# 如果接触到感染者，计算是否被感染
if n_infected_contacts > 0:
# 计算该易感者的感染概率
inf_prob = 1 - (1 - self.q) ** n_infected_contacts
# 判断是否转变为潜伏者
if np.random.random() < inf_prob:
self.stage[s_idx] = 1
new_exposed_count += 1
3.2 患病状态之间的转化
个体在潜伏期和感染期的状态转化依据持续时间进行：
1.潜伏转感染：潜伏者在潜伏期结束后转为感染者。 2.感染转康复：感染者在感染期结束后转为康复者。
更新潜伏者和感染者的状态持续时间
# 使用掩码选择潜伏者和感染者
active_mask = (self.stage == 1) | (self.stage == 2)
# 增加状态持续时间
self.duration[active_mask] += self.ts
更新潜伏者状态为感染者（潜伏期满）
# 找出潜伏期已满的潜伏者
exposed_to_infected = (self.stage == 1) & (self.duration > 1/self.sigma)
# 将其状态更新为感染者
self.stage[exposed_to_infected] = 2
# 重置其状态持续时间
self.duration[exposed_to_infected] = 0
更新感染者状态为康复者（感染期满）
# 找出感染期已满的感染者
infected_to_recovered = (self.stage == 2) & (self.duration > 1/self.gamma)
# 将其状态更新为康复者
self.stage[infected_to_recovered] = 3
4 伪代码与流程图
4.1 SEIR随机模型伪代码
| 算法1：SEIR随机模型的模拟 | |
| 输入： 总人数N，模拟时长M，时间步长ts，日接触率c，传染概率q， 潜伏期转化率sigma，感染期恢复率gamma，初始潜伏者人数E0 | |
| 输出: 每日各状态人数S_record, E_record, I_record, R_record， 每日新增潜伏者new_exposed_record，实时再生数Rt_record | |
|
1. 初始化所有个体状态(stage=0表示S，1表示E，2表示I，3表示R)； 2. 初始化每个个体在当前状态的持续时间duration； 3. for t = 1 to M do 4. 记录当前康复者人数 prev_recovered； 5. 获取当前易感者索引susceptible_indices； 6. 获取当前感染者索引infected_indices； 7. 执行易感者到潜伏者的状态转换（算法1-1）； 8. 执行潜伏者到感染者的状态转换（算法1-2）； 9. 执行感染者到康复者的状态转换（算法1-3）； 10. 记录t时刻的S, E, I, R人数； 11. 计算并记录t时刻的Rt值； 12. 记录t时刻的新增潜伏者和新增康复者人数； 13. end |
| 算法1-1：易感者到潜伏者的状态转换 | |
| 输入： 易感者索引集合susceptible_indices， 感染者索引集合infected_indices | |
| 输出: 更新后的个体状态stage | |
|
1. 初始化新增潜伏者计数器 new_exposed_count = 0； 2. if infected_indices不为空 then 3. foreach s in susceptible_indices do 4. 生成接触人数 n_contacts ~ Poisson(c * ts)； 5. if n_contacts > 0 then 6. 随机选择n_contacts个接触者contacts； 7. 计算接触感染者数量n_infected_contacts； 8. if n_infected_contacts > 0 then 9. 计算感染概率 inf_prob = 1-(1-q)^n_infected_contacts； 10. if Uniform(0,1) < inf_prob then 11. 将个体s状态更新为潜伏者(stage[s] = 1)； 12. new_exposed_count += 1； 13. end 14. end 15. end 16. end 17. end |
| 算法1-2：潜伏者到感染者的状态转换 | |
| 输入： 个体状态stage，状态持续时间duration，潜伏期转化率sigma | |
| 输出: 更新后的个体状态stage和持续时间duration | |
|
1. foreach 个体x do 2. if x.stage = 1 and x.duration > 1/sigma then 3. 将个体x状态更新为感染者(x.stage = 2)； 4. 重置个体x的状态持续时间(x.duration = 0)； 5. end 6. end |
| 算法1-3：感染者到康复者的状态转换 | |
| 输入： 个体状态stage，状态持续时间duration，感染期恢复率gamma | |
| 输出: 更新后的个体状态stage | |
|
1. foreach 个体x do 2. if x.stage = 2 and x.duration > 1/gamma then 3. 将个体x状态更新为康复者(x.stage = 3)； 4. end 5. end |
4.2 流程图
5 模拟设置
模型参数的设置对仿真结果有重要影响。以下为主要参数的设置及其含义：
-
模拟总时长（m）：100天
模拟疾病传播的总时间长度，以天为单位。 -
总人口数量（n）：10,000人
模拟的人群总数，代表了模型中的个体数量。 -
时间步长（ts）：1天
模拟的时间步长，每次迭代代表一天。 -
日接触率（c）：10人/天
每个个体每天平均接触的其他个体数量。 -
传染概率（q）：0.05
每次接触导致感染的概率。 -
潜伏期转化率（σ）：1/3天⁻¹
潜伏期转化为感染期的速率，平均潜伏期为3天。 -
感染期恢复率（γ）：1/5天⁻¹
感染期恢复为康复期的速率，平均感染期为5天。 -
初始潜伏者数量（E₀）：10人
模拟开始时，初始处于潜伏状态的个体数量。
此外，为了分析模型的不确定性，进行了多次独立模拟（如20次），并记录各次模拟的结果以进行统计分析。随机种子设置为42，以确保结果的可重复性。
if __name__ == "__main__":
# 设置模拟参数
M = 100 # 模拟时长（天）
N_SIMS = 20 # 多次模拟次数（用于SEIR和Rt的箱线图）
# 多次模拟用于SEIR和Rt的箱线图
np.random.seed(42) # 设置随机种子以保证结果可重复
models = []
for _ in range(N_SIMS):
model = SEIRRandomModel(m=M)
model.initialize_states()
model.simulate()
models.append(model)
# 绘制SEIR各类别的箱线图并连接均数折线
plot_seir_boxplot(models, M)
# 绘制Rt的箱线图并连接均数折线
plot_rt_boxplot(models, M)
# 单次模拟用于展示新增病例的随机波动
model = SEIRRandomModel(m=M)
model.initialize_states()
model.simulate()
# 绘制单次模拟的新增病例变化曲线
plot_new_cases(model)
6 模拟结果与解释
6.1模拟结果展示
通过多次模拟，生成了易感者、潜伏者、感染者、康复者数量的箱线图以及实时再生数Rt的箱线图。同时，展示了单次模拟中每日新增潜伏者和康复者的变化曲线。
1. SEIR各类别的箱线图
图1 展示了20次模拟中，易感者（S）、潜伏者（E）、感染者（I）、康复者（R）在100天内的变化分布。箱线图反映了不同时间点各类别人数的中位数、四分位数及异常值，均数则通过折线连接表示。易感者(S)数量呈单调下降趋势，从初始的10000人逐渐减少，潜伏者(E)数量在第40天左右达到峰值约1500人，随后减少，感染者(I)数量在第45天左右达到峰值约2800人，之后逐渐下降，康复者(R)数量呈S形增长，最终稳定在约9500人。
2. 实时再生数Rt的箱线图
图2 展示了20次模拟中，实时再生数Rt在100天内的变化分布。Rt的均数折线反映了疾病传播趋势的整体变化。模拟结果中，我们得到基本再生数R0 = 2.5，说明每个感染者在完全易感人群中能够将疾病传播给2.5个其他个体。随着传播过程进行，易感者数量的减少导致实时再生数Rt从初始的2.5逐渐下降。在第40天左右，Rt降至1以下，并在第60天后趋于稳定，维持在约0.2的水平，表明疫情传播强度显著减弱，但仍可能存在少量的传播。这个结果反映了在当前模型参数设置下，疾病可能会持续存在于人群中，但传播强度已经很低。
3. 单次模拟的新增病例变化曲线
图3 展示了单次模拟中，每日新增潜伏者和新增康复者的数量变化。该曲线反映了疫情发展的随机波动特征。
结果解释
-
易感者（S）：随着时间推移，易感者数量逐渐减少，特别是在疫情初期传播较快的阶段。随着更多个体转变为潜伏者和感染者，易感者的总数不断下降。
-
潜伏者（E）：潜伏者数量在疫情传播的早期快速上升，随后趋于稳定或缓慢下降。这反映了潜伏期个体的转化过程和新感染的动态平衡。
-
感染者（I）：感染者数量在潜伏者数量达到峰值后迅速上升，随后因康复者增多而逐渐下降。这表明疫情传播在达到高峰后开始受到控制。
-
康复者（R）：康复者数量持续上升，最终趋于稳定，反映了疫情结束后的免疫个体积累。
-
实时再生数Rt：
Rt的变化反映了疫情的传播力。Rt大于1表示疫情在扩散，Rt小于1表示疫情在控制。通过箱线图可以观察到Rt在不同时间点的分布及其均值的变化趋势。
6.2参数分析和讨论
总人数（𝑁）是指模型中所有个体的总数，包括易感者、潜伏者、感染者和康复者。它通常用于表示目标人群的规模。在大多数流行病模拟中，总人数决定了疾病传播的基础条件，如人群的密度、接触网络等。大规模人群流行病传播更为缓慢，峰值更高，疫情持续时间较长，群体免疫的到达需要更多的感染；小规模人群传播速度更快，峰值较高但衰退也较快。流行病通常较早达到顶峰并迅速结束。
接触频率（c）指每个个体每天与他人接触的平均次数。在此模型中，接触频率决定了易感者和感染者之间发生接触并传播的概率。较高的接触频率（如人群密集的城市环境）将大大增加病毒传播的机会。此时，尽管其他参数保持不变，病毒的传播速度会显著加快，感染者数量的增长也会更为迅速；在接触频率较低的情况下（如社交距离较远或密闭环境中），疾病的传播会变慢，整个流行病的峰值也会降低。
传染期（invGamma）是指个体从感染状态转为康复状态所需的时间，通常以天为单位，传染期的倒数通常称为康复率（γ）。如果传染期较长，意味着感染者在传染阶段停留的时间较长，这将增加每个感染者传播病毒的机会，从而可能导致更高的感染人数；较短的传染期则意味着感染者快速康复，减少了传播的机会，可能会导致较低的感染人数，尤其在疾病传播速度较快时，整个流行病会较早达到峰值并迅速消退。
平均潜伏期（invSigma）是指个体从暴露于病原体到出现症状并具有传染性所需的时间，通常以天为单位，潜伏期的倒数通常称为暴露率（σ）。当潜伏期较短时，疾病在人群中传播迅速，易感者数量快速减少，潜伏者和感染者人数会在短期内达到较高的峰值。而当潜伏期延长时，疾病传播速度明显放缓，易感者数量下降更加缓慢，潜伏者和感染者的峰值显著降低但分布更加平缓，康复者数量增长缓慢，使得整个流行过程明显延长。
7 源代码
模拟用时9.57秒
import numpy as np
import matplotlib.pyplot as plt
class SEIRRandomModel:
def __init__(self, m=100, n=10000, ts=1, c=10, q=0.05, sigma=1/3, gamma=1/5):
"""
初始化SEIR随机模型参数
Args:
m (int): 模拟总时长（天）
n (int): 总人口数量（人）
ts (float): 时间步长（天）
c (float): 日接触率（人/天）
q (float): 传染概率（无量纲）
sigma (float): 潜伏期转化率（1/天）
gamma (float): 感染期恢复率（1/天）
"""
self.m = m # 模拟总时长（天）
self.n = n # 总人口数量（人）
self.ts = ts # 时间步长（天）
self.c = c # 日接触率（人/天）
self.q = q # 传染概率（无量纲）
self.sigma = sigma # 潜伏期转化率（1/天）
self.gamma = gamma # 感染期恢复率（1/天）
def initialize_states(self, initial_exposed=10):
"""
初始化模型状态
Args:
initial_exposed (int): 初始潜伏者数量（人）
"""
# 初始化个体状态数组（0:S, 1:E, 2:I, 3:R）
self.stage = np.zeros(self.n, dtype=int)
self.stage[0:initial_exposed] = 1 # 设置初始潜伏者
self.duration = np.zeros(self.n) # 记录个体在当前状态的持续时间（天）
# 初始化时间序列记录数组
self.S_record = np.zeros(self.m) # 易感者数量记录（人）
self.E_record = np.zeros(self.m) # 潜伏者数量记录（人）
self.I_record = np.zeros(self.m) # 感染者数量记录（人）
self.R_record = np.zeros(self.m) # 康复者数量记录（人）
self.Rt_record = np.zeros(self.m) # 实时再生数记录（无量纲）
# 新增记录数组
self.new_exposed_record = np.zeros(self.m) # 新增潜伏者数量记录（人/天）
self.new_recovered_record = np.zeros(self.m) # 新增康复者数量记录（人/天）
def simulate(self):
"""运行SEIR模型仿真"""
for i in range(self.m):
# 记录本时间步开始时的状态，用于计算新增人数
prev_recovered = np.sum(self.stage == 3)
# 获取当前易感者和感染者的索引
susceptible_indices = np.where(self.stage == 0)[0]
infected_indices = np.where(self.stage == 2)[0]
# 记录新增潜伏者的计数器
new_exposed_count = 0
if len(infected_indices) > 0: # 只在有感染者时进行传播计算
# 为每个易感者生成接触者
for s_idx in susceptible_indices:
# 生成该易感者当天接触的人数（泊松分布）
n_contacts = np.random.poisson(self.c * self.ts)
if n_contacts > 0:
# 随机选择接触者
contacts = np.random.choice(self.n, size=n_contacts)
# 计算接触的感染者数量
n_infected_contacts = np.sum(np.isin(contacts, infected_indices))
# 如果接触到感染者，计算是否被感染
if n_infected_contacts > 0:
# 计算该易感者的感染概率
inf_prob = 1 - (1 - self.q) ** n_infected_contacts
# 判断是否转变为潜伏者
if np.random.random() < inf_prob:
self.stage[s_idx] = 1
new_exposed_count += 1
# 更新潜伏者和感染者的状态持续时间
active_mask = (self.stage == 1) | (self.stage == 2)
self.duration[active_mask] += self.ts
# 更新潜伏者状态为感染者（潜伏期满）
exposed_to_infected = (self.stage == 1) & (self.duration > 1/self.sigma)
self.stage[exposed_to_infected] = 2
self.duration[exposed_to_infected] = 0
# 更新感染者状态为康复者（感染期满）
infected_to_recovered = (self.stage == 2) & (self.duration > 1/self.gamma)
self.stage[infected_to_recovered] = 3
# 记录各类人群数量和再生数
self.S_record[i] = np.sum(self.stage == 0)
self.E_record[i] = np.sum(self.stage == 1)
self.I_record[i] = np.sum(self.stage == 2)
self.R_record[i] = np.sum(self.stage == 3)
self.Rt_record[i] = (self.c * self.q) / self.gamma * (self.S_record[i] / self.n)
# 记录新增潜伏者和康复者数量
self.new_exposed_record[i] = new_exposed_count
self.new_recovered_record[i] = np.sum(self.stage == 3) - prev_recovered
def plot_seir_boxplot(models, m):
"""
绘制SEIR模型多次模拟的每日箱线图，并连接每日均数的折线
Args:
models (list): SEIR模型对象列表
m (int): 模拟总时长（天）
"""
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 创建四个子图
fig, axes = plt.subplots(2, 2, figsize=(20, 15))
titles = ['易感者(S)', '潜伏者(E)', '感染者(I)', '康复者(R)']
records = ['S_record', 'E_record', 'I_record', 'R_record']
colors = ['g', 'y', 'r', 'b']
t = np.arange(0, m, 1) # 每一天画一个箱线图
for idx, (title, record, color, ax) in enumerate(zip(titles, records, colors, axes.flat)):
# 收集每个时间点的数据
data = [ [getattr(model, record)[time] for model in models] for time in t]
# 绘制箱线图
ax.boxplot(data, positions=t, widths=0.6, patch_artist=True,
boxprops=dict(facecolor=color, alpha=0.5))
ax.set_title(title, fontsize=16)
ax.set_xlabel('时间（天）', fontsize=14)
ax.set_ylabel('人数（人）', fontsize=14)
ax.grid(True, linestyle='--', alpha=0.7)
# 设置x轴刻度，每隔五天一个刻度，并增大字体
ax.set_xticks(np.arange(0, m, 5))
ax.set_xticklabels(np.arange(0, m, 5), rotation=45, fontsize=12)
# 计算每日均数
means = [np.mean([getattr(model, record)[time] for model in models]) for time in t]
# 绘制均数折线
ax.plot(t, means, color='k', linewidth=2, label='均数')
ax.legend()
plt.tight_layout()
plt.show()
def plot_rt_boxplot(models, m):
"""
绘制Rt多次模拟的每日箱线图，并连接每日均数的折线
Args:
models (list): SEIR模型对象列表
m (int): 模拟总时长（天）
"""
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.figure(figsize=(20, 10))
t = np.arange(0, m, 1) # 每天画一个箱线图
# 收集每个时间点的Rt数据
data = [ [model.Rt_record[time] for model in models] for time in t ]
# 绘制箱线图
plt.boxplot(data, positions=t, widths=0.6, patch_artist=True,
boxprops=dict(facecolor='cyan', alpha=0.5))
# 设置x轴刻度，每隔五天一个刻度，并增大字体
plt.xticks(np.arange(0, m, 5), rotation=45, fontsize=12)
# 计算Rt的每日均数
rt_means = [np.mean([model.Rt_record[time] for model in models]) for time in t]
# 绘制Rt均数折线
plt.plot(t, rt_means, color='k', linewidth=2, label='Rt均数')
plt.title('实时再生数Rt变化分布', fontsize=18)
plt.xlabel('时间（天）', fontsize=16)
plt.ylabel('Rt（无量纲）', fontsize=16)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=14)
plt.tight_layout()
plt.show()
def plot_new_cases(model):
"""
绘制单次模拟的新增潜伏者和新增康复者曲线
Args:
model: 单个SEIR模型对象
"""
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.figure(figsize=(12, 8))
t = np.arange(model.m)
# 直接绘制单次模拟的结果
plt.plot(t, model.new_exposed_record, 'r-', label='新增潜伏者', linewidth=2)
plt.plot(t, model.new_recovered_record, 'g-', label='新增康复者', linewidth=2)
plt.title('每日新增病例数变化', fontsize=16)
plt.xlabel('时间（天）', fontsize=14)
plt.ylabel('人数（人/天）', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=12, loc='upper right')
plt.xticks(np.arange(0, model.m, 10), fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
if __name__ == "__main__":
# 设置模拟参数
M = 100 # 模拟时长（天）
N_SIMS = 20 # 多次模拟次数（用于SEIR和Rt的箱线图）
# 多次模拟用于SEIR和Rt的箱线图
np.random.seed(42) # 设置随机种子以保证结果可重复
models = []
for _ in range(N_SIMS):
model = SEIRRandomModel(m=M)
model.initialize_states()
model.simulate()
models.append(model)
# 绘制SEIR各类别的箱线图并连接均数折线
plot_seir_boxplot(models, M)
# 绘制Rt的箱线图并连接均数折线
plot_rt_boxplot(models, M)
# 单次模拟用于展示新增病例的随机波动
model = SEIRRandomModel(m=M)
model.initialize_states()
model.simulate()
# 绘制单次模拟的新增病例变化曲线
plot_new_cases(model)