SEIR 模型的参数如何确立，接触率动态变化

Question

https://kns.cnki.net/KXReader/Detail?invoice=Kv7yhz1FnWjxLqDjBGhtK1ECGs1VQ%2FQPNJUBcTgOd9O%2F9DGN%2BKjn29QLRsxl3b5L1o56bNUJMnr8m3LIT259%2BJvZ1fSpFmQaA8GRITb3RezQff%2FFtaADv07zQjcDi%2Bq7KKvx06vWSvaUPQgosbe2kHJ5P67zGWL078Paalsbk9k%3D&DBCODE=CJFQ&FileName=XDXK202323030&TABLEName=cjfdlast2024&nonce=C9D9AFC99E454DE193EDEA3599B73F9D&TIMESTAMP=1713923538243&uid=
基于动态 SEIR 模型的小规模疫情预测 ，这篇文章，里面的这个式子是怎么得出来的呢？

Answer 1

引自免费微信小程序：皆我百晓生

在理解动态SEIR模型的参数确立过程中，我们首先要了解传统SEIR模型的基本构成。SEIR模型是一个流行病学模型，描述了人群中的四种状态：易感者(Susceptible)，暴露者(Exposed)，感染者(Infectious)，和康复者(REmoved)。在传统模型中，感染率(β)和恢复率(γ)通常被视为常数。

然而，在实际应用中，特别是在 COVID-19 这样的传染病中，接触率（即感染率）可能随着时间和外部因素（例如防控措施、季节性变化等）而改变。文章提到的方法是引入一个动态的接触率模型，以便更好地反映实际情况。

文中提出的动态接触率模型使用的是逻辑斯蒂函数 (Logistic Function)，这是一种常用的数学模型，用来描述某个量在有限时间内从低值向高值增长的过程，其形状类似于"S"形曲线。逻辑斯蒂函数可以很好地模拟人口行为的变化，比如防疫意识的提升或下降。

具体来说，动态接触率 Ce 可以表示为：

[ C_e = \frac{C}{\left(1 + e^{-a(t+m)}\right)^2} ]

这里，

• ( C ) 是一个尺度常数，决定了最大接触率的水平。
• ( a ) 是正则化参数，影响曲线的陡峭程度，也就是反应速度。
• ( t ) 是时间变量，从发现首例感染病例开始计时。
• ( m ) 是偏置参数，调整曲线的位置，可以对应于防控措施实施的时间点。

为了确定这些参数的最优值，文章采用了机器学习中的优化方法，特别是最小化均方误差(Mean Squared Error, MSE)作为损失函数。MSE 量化了模型预测值与真实数据之间的差异平方的平均值，目标是使这个差异尽可能小。对于给定的数据集，我们会找到一组参数 ( a, C, m )，使得预测的接触率与实际观察到的接触率之间的MSE最小。

具体的计算步骤包括以下部分：

1. 数据准备：收集关于疫情发展的数据，包括每个时间点的确诊病例数量或其他能间接反映接触率的数据。
2. 设定初始参数：对 ( a, C, m ) 设置初步猜测值。
3. 拟合模型：使用这些参数计算每一天的动态接触率，并与实际数据比较。
4. 损失函数：计算损失函数，即预测值与实际值之间差的平方和。
5. 优化参数：使用梯度下降法或者其他优化算法，调整 ( a, C, m ) 来最小化损失函数。
6. 迭代：重复步骤3到5，直到参数收敛，即损失函数不再显著减小。
7. 评估模型：使用优化后的参数重新运行模型，并评估其预测性能。

这是一段简化的Python代码示例，展示如何使用scipy库中的optimize.minimize函数来寻找最佳参数：

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# 假设已有的数据
true_contact_rates = [your_data_here]  # 时间序列的真实接触率数据
initial_infections_time = your_initial_infection_time  # 发现首例感染病例的时间点

def logistic_function(t, a, C, m):
    return C / ((1 + np.exp(-a * (t - m)))**2)

def loss_function(params, true_contact_rates, initial_infections_time):
    a, C, m = params
    predicted_contact_rates = [logistic_function(i, a, C, m) for i in range(initial_infections_time, len(true_contact_rates))]
    return sum((predicted_contact_rates - true_contact_rates)**2)

# 初始参数
params = np.array([0.1, 10, 0])

# 使用最小化损失函数的算法
solution = minimize(loss_function, params, args=(true_contact_rates, initial_infections_time))

# 输出优化后的参数
optimized_params = solution.x
print("Optimized parameters: a =", optimized_params[0], ", C =", optimized_params[1], ", m =", optimized_params[2])

# 使用优化后的参数绘制预测曲线
predicted_contact_rates = [logistic_function(i, *optimized_params) for i in range(initial_infections_time, len(true_contact_rates))]
plt.plot(range(len(true_contact_rates)), true_contact_rates, label="True Contact Rates")
plt.plot(range(initial_infections_time, len(true_contact_rates)), predicted_contact_rates, label="Predicted Contact Rates")
plt.xlabel("Time")
plt.ylabel("Contact Rate")
plt.legend()
plt.show()

请注意，这段代码仅作演示用途，你需要根据实际情况替换数据和参数。此外，优化过程可能会受到初始化参数的影响，因此可能需要尝试不同的初始值以获得较好的结果。同时，确保数据的质量和完整性对优化过程至关重要。