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模型（Susceptible-Exposed-Infectious-Recovered）是一种流行病学模型，用来描述传染病在人群中的传播过程。在动态SEIR模型中，接触率（contact rate）确实会随着时间和政策干预而变化。

公式 Ce - a(t + m) / (1 + e^-a(t + m))^2 来自于逻辑斯蒂函数（sigmoid function）的变形，通常用于表示接触率随时间的非线性增长和衰减。这里的参数解释如下：

• Ce 是接触率的最大值（饱和值），也就是在无任何限制时可能达到的最高接触水平。
• a 是正则化参数，它决定了接触率变化的速度。较大的 a 值意味着接触率变化更快。
• t 是发现首例感染的时间点。
• m 是偏置，表示接触率从 Ce 开始变化的滞后时间。

公式中的 (1 + e^-a(t + m))^2 是逻辑斯蒂函数的倒数，通常用于标准化，确保接触率始终在0到1之间，即使 t 很大，接触率也不会无限增长。

均方损失函数（Mean Squared Error, MSE）被用来评估模型预测的准确性，即实际接触率与模型预测接触率的差异的平方的平均值。通过最小化MSE，可以找到最优的 a、C 和 m 参数值。

完整的参数估计过程可以通过优化算法实现，比如梯度下降法或拟牛顿法，结合Python的Scipy库中的optimize.minimize函数。这里是一个简化版的代码示例，假设你有一个包含实际接触率数据的时间序列 observed_contacts：

import numpy as np
from scipy.optimize import minimize

def objective_function(params, t, observed_contacts):
    a, C, m = params
    predicted_contacts = C * (1 + np.exp(-a * (t + m))) ** (-2)
    return np.mean((predicted_contacts - observed_contacts) ** 2)

# 初始化参数猜测
initial_params = [0.5, 100, 10]  # a, C, m
t_data = np.arange(start=0, stop=len(observed_contacts), step=1)  # 时间序列
observed_contacts = ...  # 你的实际接触率数据

# 使用Scipy的minimize函数进行优化
optimization_result = minimize(objective_function, initial_params, args=(t_data, observed_contacts))

# 获取最优参数
optimal_a, optimal_C, optimal_m = optimization_result.x

# 使用最优参数更新模型
dynamic_contact_rate = [optimal_C * (1 + np.exp(-optimal_a * (t + optimal_m))) ** (-2) for t in t_data]

这个代码只是一个基本框架，实际应用中你可能需要根据具体的数据格式和实际情况调整。记得在运行之前确保你有实际的接触率数据 observed_contacts，并根据需要对数据进行预处理。如果你有实际数据，我可以帮你进一步细化代码。