ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

Python疫情数据分析源码复现:SEIR模型与动态接触率实战

Python疫情数据分析源码复现:SEIR模型与动态接触率实战 简介这份资源是面向数据分析初学者、科研人员与公共卫生工作者的COVID-19疫情数据分析与可视化源码合集基于Python与Jupyter Notebook实现可用于研究疫情传播模式、发展趋势与高风险区域也适合作为数据科学入门练手项目。压缩包共66个文件约8.6MB包含15个ipynb交互式分析笔记、15张png可视化结果图、12个py脚本、6个csv疫情数据文件以及pkl模型文件、txt说明与pdf文档等覆盖数据处理、建模预测与图表展示的完整流程。目录按Model 1、Model 2、Model 3等模块划分涉及SEIR模型、动态SEIR预测、岭回归与多项式回归基线等具体分析主题结构清晰便于按模块查阅。目前已有333人学习下载读者可借此掌握从数据清洗、模型构建到结果可视化的完整思路并直接复用脚本与笔记快速开展自己的疫情分析实践。1. 从一份 70 文件的源码包说起COVID-19 数据分析到底能复现出什么2020 年那会儿很多人第一次意识到「数据能决定一件事的走向」。这份基于 Python 的 COVID-19 疫情数据分析和可视化源码包就是那个阶段留下来的一套完整工程70 个文件15 个 Jupyter Notebook、12 个 Python 脚本、6 个 CSV 数据文件、15 张 PNG 结果图外加 requirements.txt 和 readme.txt。它不是一篇论文的附件而是一个能跑起来、能改参数、能换数据源的分析工作台。核心围绕三块Model 1 用报告数据反推武汉真实感染规模Model 2 用 SEIR/SIR 做疫情趋势预测Model 3 用动态 SEIR 加 Ridge 回归、多项式回归做基线对比。适合公共卫生方向的研究生、做数据分析项目练手的工程师以及需要一套可复用疫情建模模板的从业者。下面我按「能跑通 → 能改 → 不翻车」的顺序拆一遍。2. 环境与数据链路把 Notebook 从「打不开」跑到「出图」2.1 依赖安装与 Jupyter 启动的最小闭环拿到压缩包先别急着双击 ipynb。这个项目里有 requirements.txt但疫情类项目常见的坑是当年写的 pandas、numpy、matplotlib 版本和现在差了好几个大版本直接pip install -r大概率会在某个pd.append或者sns.distplot上报错。我一般会先建独立虚拟环境再按需降级。# 建环境Python 3.7 附近最稳因为源码里有 cpython-37 的 pyc python -m venv venv_covid source venv_covid/bin/activate # Windows 用 venv_covid\Scripts\activate # 先装核心三件套版本别追新 pip install numpy1.21.6 pandas1.3.5 matplotlib3.5.3 pip install jupyter notebook scikit-learn1.0.2 scipy1.7.3 # 再补项目可能用到的 pip install seaborn plotly statsmodels逻辑说明源码目录里出现了metadata.cpython-37.pyc说明作者当年跑在 Python 3.7用 3.7~3.9 兼容性最好。参数上pandas 卡在 1.3.x 是因为 2.0 之后删掉了append、改了inplace行为很多老 Notebook 会直接崩。装完执行jupyter notebook浏览器打开后先别跑全部 cell逐个 Kernel → Restart Run All 之前先确认数据路径。2.2 数据文件与路径DXY_Chinese.csv 和 DXYArea.csv 怎么接项目 data 目录下有两个关键 CSVDXY_Chinese.csv和DXYArea.csv这是早期丁香园DXY口径的疫情数据。Notebook 里读数据的写法通常是相对路径一旦你把 ipynb 挪了位置就会FileNotFoundError。稳妥做法是在每个 Notebook 开头统一加一段路径锚定。import os import pandas as pd # 以项目根目录为基准避免相对路径漂移 BASE_DIR os.path.dirname(os.path.abspath(__file__)) DATA_DIR os.path.join(BASE_DIR, data) df_area pd.read_csv(os.path.join(DATA_DIR, DXYArea.csv)) df_cn pd.read_csv(os.path.join(DATA_DIR, DXY_Chinese.csv)) # 先看列名和日期范围别急着画图 print(df_area.columns.tolist()) print(df_area[updateTime].min(), df_area[updateTime].max())逻辑说明os.path.abspath(__file__)在 Notebook 里拿不到真实文件路径更稳的是用os.getcwd()并约定「从项目根目录启动 jupyter」。参数上重点看updateTime字段DXY 数据是累积快照同一天可能有多条画趋势前必须按provinceName、updateTime去重取最新。DXY_Chinese.csv偏中文省市维度DXYArea.csv偏地理区域维度两张表 join 之前先统一省市英文名——项目里那两个 pklchineseProvince_to_EN.pkl、chineseCity_to_EN.pkl就是干这个的。2.3 三个 Model 的职责边界Model 1 的Estimating_current_cases_in_Wuhan.ipynb做的是「从境外输入病例反推武汉当时真实感染数」属于回溯估计Model 2 的Forecast_Outbreak_Wuhan.ipynb用 SEIR/SIR 做未来趋势Model 3 的Forecast_by_DynamicSEIR.ipynb把接触率做成随时间变化的动态参数再和 Ridge、多项式回归做基线对比。三者不是递进关系而是三种独立视角别指望跑完 Model 1 自动喂给 Model 2数据要自己对齐。3. SEIR 与动态接触率参数怎么设、图怎么读3.1 SEIR_model.py 的四个仓室与关键参数SEIR 把人群分成易感S、暴露E、感染I、移除R四类。项目里SEIR_model.py和helper_fun_epi_model.py是核心前者定义微分方程后者提供拟合和绘图辅助。常见做法是用scipy.integrate.odeint求解。import numpy as np from scipy.integrate import odeint def seir_deriv(y, t, beta, sigma, gamma, N): S, E, I, R y # beta 接触率sigma 潜伏转感染率gamma 移除率 dS -beta * S * I / N dE beta * S * I / N - sigma * E dI sigma * E - gamma * I dR gamma * I return dS, dE, dI, dR N 11_000_000 # 武汉常住人口量级 beta, sigma, gamma 0.6, 1/5.2, 1/7.0 y0 [N-1, 0, 1, 0] t np.linspace(0, 120, 120) sol odeint(seir_deriv, y0, t, args(beta, sigma, gamma, N))逻辑说明beta是单位时间有效接触率直接决定峰值高度sigma取 1/潜伏期常见潜伏期 5.2 天gamma取 1/感染期约 7 天。这三个参数不是拍脑袋项目里通过拟合早期确诊曲线反推。N用武汉人口量级换城市必须改否则 S 的下降速度会失真。跑完sol四列分别对应 S/E/I/R画I列就是感染人数曲线。3.2 动态 SEIR把 beta 变成时间的函数Model 3 的Dynamic_SEIR_model.py比静态版多做了一件事让接触率beta随干预措施变化。项目里contact_rate.png、beta.png两张图就是它的输出。实现上通常把 beta 写成阶跃或指数衰减函数。def beta_t(t, beta0, t_lock, decay): # t_lock 之前保持 beta0之后按 decay 衰减模拟封控效果 if t t_lock: return beta0 return beta0 * np.exp(-decay * (t - t_lock)) # 在 odeint 的 deriv 里按当前 t 取 beta def seir_dynamic(y, t, beta0, t_lock, decay, sigma, gamma, N): beta beta_t(t, beta0, t_lock, decay) S, E, I, R y dS -beta * S * I / N dE beta * S * I / N - sigma * E dI sigma * E - gamma * I dR gamma * I return dS, dE, dI, dR逻辑说明t_lock是干预生效的时间点decay越大代表措施越狠、接触率掉得越快。这两个参数是动态 SEIR 的灵魂调它们能明显改变峰值出现的时间和高度。项目里SEIR_prediction.png、SEIR_test_7days.png就是拿预测值和真实值对比看模型有没有跑偏。注意odeint要求 deriv 函数签名固定把额外参数用args传进去别用全局变量否则换参数要改代码。3.3 基线模型Ridge 与多项式回归为什么必须留着Model 3 里Baseline_RidgeRegression.ipynb和Baseline_polynomial_regression.ipynb不是凑数。疫情预测里复杂模型很容易过拟合一个简单的 Ridge 回归反而能当「下限参照」。项目输出的baseline_ridge.png、residuals_ridge_train_0216.csv、residuals_ridge_test_0216.csv就是拿残差说话。from sklearn.linear_model import Ridge from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import make_pipeline # 用前 N 天特征预测后一天新增 X_train, y_train build_lag_features(df_cn, lag7) model make_pipeline(PolynomialFeatures(2), Ridge(alpha1.0)) model.fit(X_train, y_train) # 看训练/测试残差判断有没有系统性偏差 residuals y_train - model.predict(X_train) print(residuals.describe())逻辑说明lag7表示用过去 7 天数据构造特征alpha1.0是 Ridge 的正则强度越大越保守。残差 CSV 是项目已经存好的直接读出来画residuals_ridge_test_0216.csv的分布如果残差均值明显偏离 0说明模型有系统偏差动态 SEIR 的结果也要打个问号。这一步是很多人会跳过的验证环节但恰恰是判断「模型能不能信」的关键。4. 可视化与数据清洗15 张 PNG 背后的坑4.1 中文省市名映射与去重项目里chineseCity_to_EN.pkl、chineseProvince_to_EN.pkl两个 pickle 是清洗利器。DXY 数据里省市名有中文、有简写、有「境外输入」这种非地理条目直接 groupby 会得到一堆脏分组。import pickle with open(data_processing/chineseProvince_to_EN.pkl, rb) as f: prov_map pickle.load(f) df_area[province_EN] df_area[provinceName].map(prov_map) # 映射不上的先看有哪些别直接 drop unmapped df_area[df_area[province_EN].isna()][provinceName].unique() print(unmapped)逻辑说明map之后一定要检查isna()因为「境外输入」「待明确地区」这类条目本来就不该进省级趋势图。常见做法是先打印未映射清单人工确认后再决定是丢弃还是单独归类。去重时按province_ENupdateTime取最后一条因为 DXY 是累积快照同一天多次更新只保留最新。4.2 累积值转新增值别在累积曲线上做拟合DXYArea.csv里的confirmedCount是累积确诊很多 Notebook 直接拿它画趋势、做回归这是典型误用。累积值单调递增任何回归都会得到「一直涨」的假象。正确做法是先差分。df_area df_area.sort_values([province_EN, updateTime]) df_area[new_confirmed] ( df_area.groupby(province_EN)[confirmedCount].diff().fillna(0) ) # 差分可能出现负数数据修正置零处理 df_area.loc[df_area[new_confirmed] 0, new_confirmed] 0逻辑说明groupby().diff()按省做一阶差分得到每日新增fillna(0)处理每个省的第一天。负数来自数据回溯修正置零是常见妥协但要在报告里注明。项目里Incorrect_confirmed_cases.ipynb专门处理这类异常值得先跑一遍再动主分析。4.3 出图规范15 张 PNG 的命名与复现项目 image 目录下withControl.png、without_control.png、dynamic_SEIR.png、SEIR_fit.png等图命名已经说明了对比关系。复现时建议固定 matplotlib 的dpi和figsize否则出的图和原图对不上评审时会被质疑。import matplotlib.pyplot as plt plt.rcParams[figure.dpi] 120 plt.rcParams[savefig.dpi] 120 plt.rcParams[font.size] 11 fig, ax plt.subplots(figsize(9, 5)) ax.plot(t, sol[:, 2], labelInfected (I)) ax.plot(t, sol[:, 0], labelSusceptible (S)) ax.set_xlabel(Days); ax.set_ylabel(Population) ax.legend(); fig.tight_layout() fig.savefig(image/seir_reproduced.png)逻辑说明dpi120是折中值太低糊、太高文件大。font.size统一避免图注大小不一。保存路径对齐项目 image 目录方便和原图并排比对。如果中文字体报错加plt.rcParams[axes.unicode_minus] False并指定系统中文字体。5. 避坑与排查跑这份源码最容易翻车的五处现象一Notebook 打开后所有 cell 报ModuleNotFoundError: No module named helper_fun_epi_model。原因Notebook 和 helper 脚本不在同一目录或者启动 jupyter 时的工作目录不是项目根目录。 解决cd到项目根目录再jupyter notebook或在 Notebook 开头加sys.path.append(os.path.abspath(..))把 Model 目录加进搜索路径。现象二SEIR 拟合出的曲线和SEIR_fit.png差很远峰值高得离谱。原因beta或N没按目标城市改或者odeint的时间步长太粗导致数值发散。 解决先确认N是目标城市人口量级再把t np.linspace(0, 120, 120)加密到 1200 个点beta从 0.3 开始逐步试别一上来就 0.6。现象三Ridge 回归残差 CSV 读出来全是 NaN。原因residuals_ridge_test_0216.csv里可能有空行或索引列错位pandas 默认把首列当索引。 解决pd.read_csv(path, index_col0)或pd.read_csv(path).dropna()先head()看结构再决定。现象四中文省市映射后大量 NaN趋势图只剩几个省。原因DXY 数据版本更新后省市名变了旧 pkl 映射表覆盖不全。 解决打印未映射清单手动补映射字典别直接dropna()否则会丢掉「境外输入」等关键类别导致总量对不上。现象五动态 SEIR 的beta_t函数在odeint里报参数数量错误。原因odeint要求 deriv 函数签名是(y, t, *args)把beta0、t_lock等塞进全局变量或闭包容易出错。 解决统一用args(beta0, t_lock, decay, sigma, gamma, N)传参函数定义里按顺序接收别混用全局变量。6. 进阶把这份源码改成自己的分析模板跑通只是第一步真正有价值的是把它变成可复用的模板。我一般会做三件事。第一把data_processing/DXY_AreaData_query.py里的数据读取逻辑抽成一个load_data()函数统一处理路径、去重、差分所有 Notebook 都调它避免每个文件重复写清洗代码。第二把 SEIR 的参数拟合封装成fit_seir(df, city, n_days)输入一个省市的每日新增序列输出最优beta、sigma、gamma和拟合图这样换城市只要改一个参数。第三给动态 SEIR 加一个参数扫描。import itertools import numpy as np results [] for beta0, decay in itertools.product([0.4, 0.5, 0.6], [0.02, 0.05, 0.1]): sol odeint(seir_dynamic, y0, t, args(beta0, 50, decay, sigma, gamma, N)) peak sol[:, 2].max() peak_day t[sol[:, 2].argmax()] results.append({beta0: beta0, decay: decay, peak: peak, peak_day: peak_day}) import pandas as pd print(pd.DataFrame(results).sort_values(peak))逻辑说明itertools.product做笛卡尔积扫描t_lock50固定为干预生效日输出每个参数组合下的峰值和峰值日。这张表能直接告诉你「接触率降得越快峰值越低、来得越晚」比单张曲线图更有说服力。参数上decay从 0.02 到 0.1 覆盖了温和到严格干预beta0从 0.4 到 0.6 覆盖常见估计区间。验证方法上我习惯留一段「后悔药」把最后 7 天数据从训练集里切出来当 holdout跑完预测再对比SEIR_test_7days.png的思路算 MAE 和 MAPE。如果 MAPE 超过 30%说明参数过拟合得回去调sigma、gamma或者换动态模型。这套流程我从那以后每次做疫情类建模都强制走一遍——先差分、再拟合、留 holdout、扫参数四步缺一不可。希望帮到你。本文还有配套的精品资源点击获取
返回列表