简介一套Python实现的SIR模型模拟项目聚焦BA无标度网络与ER随机网络上的传染病传播对比。面向网络科学、流行病学建模初学者以及Python数据分析学习者可直观理解网络结构对疾病扩散的影响。压缩包共21个文件包含4个Python脚本主程序、BA/ER网络构建与SIR模拟、10张仿真结果图、README说明文档及配置附件总大小仅217KB轻量便于阅读与二次开发。代码覆盖SIR易感-感染-康复状态迁移、网络生成与参数可视化结合输出图表可观察BA网络中高连接节点导致的集中爆发以及ER网络中相对均匀的传播特征。README与文件组织清晰适合作为网络动力学实验的起步模板。已有875人学习下载对理解真实复杂网络中的疫情传播规律具有实践参考价值。1. 为什么要在 BA/ER 网络上做 SIR 仿真一个反直觉的传播阈值结论我第一次用这个方案跑通仿真时最反直觉的结论是在 BA 无标度网络上SIR 模型的传播阈值几乎为零在 ER 随机网络上阈值清晰可见。同样低的传染率疾病能在 BA 网络上形成大流行却很难在 ER 网络上扩散。原因是度分布结构不同BA 网络靠择优连接生成存在少数度极高的 hub 节点ER 网络的度集中在平均值附近。标题里的 BA_ER_SIR_python-master 就是把这两种网络与 SIR 模型组合在一起的 Python 仿真项目回答“网络结构如何决定疫情爆发规模与速度”。适合网络科学课程作业、论文复现以及想评估真实网络传播风险的工程师不需要额外平台装好 Python、networkx、numpy 就能跑。2. 网络构建先让 BA 和 ER 两张图在平均度上对齐2.1 BA 无标度网络的生成原理择优连接与幂律度分布BA 网络的标准生成算法是 Barabási–Albert 模型从一个小规模初始图开始每次新增一个节点让它带出 m 条边这些边连接到已有节点的概率正比于该节点的当前度。度大的节点被选中的概率越来越大形成“富者愈富”的正反馈。最终生成的网络里度分布近似满足 P(k) ~ 2m² k⁻³也就是幂律指数约为 3 的无标度网络。这个幂律带来两条关键特性一是有一小撮 hub 节点的度远高于平均水平二是网络里不存在一个“典型度值”供你描述大多数节点。正因如此它才叫无标度。ER 网络的生成模型是 G(n, p)固定 n 个节点每对节点以概率 p 独立连边。度分布是二项分布网络规模够大时趋近泊松分布节点度都集中在均值附近最大的那个节点也不至于大到哪里去。把这两个网络放在一起对比本质上是比较“异质网络”与“同质网络”在传播动力学上的行为差异。这也是这个项目最核心的看点相同的平均度、相同的节点数仅因为度的方差不同传染病就会走出完全不同的曲线。选这两个网络还有一个实操上的好处它们的平均度都能直接由参数算出来。BA 网络在 n 较大时平均度近似等于 2mER 网络的平均度按定义是 p(n-1)。把两个平均度对齐之后所有传播结果差异才能归因于“度分布结构”而不是“边的多少”。2.2 用 networkx 生成网络并检查度分布log-log 图上是否成直线先装好依赖pip install networkx numpy。然后按下面的代码生成两张图。这里最关键的是 p 的计算令 ER 的期望平均度 p(N-1) 等于 BA 的平均度 2m反推得到 p 2m/(N-1)。在这个设定下BA 的边数约 2991ER 的期望边数约 3000几乎一致。如果随手把 p 设成 0.1ER 的边数会多一个量级传播速度和规模就会被“边更多”主导后面的对比结论全站不住。import numpy as np import networkx as nx N 1000 # 节点数 m 3 # BA 新增节点的连边数平均度约 2m6 p 2 * m / (N - 1) # ER 连边概率让平均度对齐 G_ba nx.barabasi_albert_graph(N, m) G_er nx.erdos_renyi_graph(N, p) print(BA 边数:, G_ba.number_of_edges()) print(ER 边数:, G_er.number_of_edges()) ba_deg np.array([d for _, d in G_ba.degree()]) er_deg np.array([d for _, d in G_er.degree()]) print(fBA 平均度 {ba_deg.mean():.3f} 最大度 {ba_deg.max()}) print(fER 平均度 {er_deg.mean():.3f} 最大度 {er_deg.max()})注意barabasi_albert_graph的第二个参数是新增节点的连边数 m不是最终平均度。m 越小网络越稀疏hub 相对越突出实验里 m 取 2 到 5 都常见。网络生成对了没有光靠边数和平均度不够还要看度分布的形状。我一般会写一个通用的度分布统计函数把每个度值的节点占比算出来后面换网络也能复用def degree_distribution(G): deg np.array([d for _, d in G.degree()]) hist np.bincount(deg).astype(float) return hist / hist.sum() hist_ba degree_distribution(G_ba) hist_er degree_distribution(G_er) print(BA 度分布前 10 项:, hist_ba[:10]) print(ER 最大度位置:, hist_er.argmax())np.bincount返回的是“每个度值的节点个数”除以总数就是概率质量函数数组下标本身就是度值读起来很直观。BA 网络的最小度是 m3因此度 0、1、2 的概率为 0属于正常现象。接着用 matplotlib 画图BA 用 log-log 坐标尾部近似一条下降直线说明幂律成立ER 则在均值附近形成一个明显的单峰鼓包不会出现长尾。两条曲线形态差异就是后面所有动力学差异的结构根源。2.3 平均度对齐的验证与邻接矩阵导出实验脚本里我习惯加一道断言防止中途改了某个参数后两张图悄悄失配。先验证平均度再把邻接矩阵转成 numpy 数组供后续 SIR 模块使用assert abs(ba_deg.mean() - er_deg.mean()) 0.2, 平均度偏差过大检查 p 的取值 A_ba nx.to_numpy_array(G_ba).astype(np.int8) A_er nx.to_numpy_array(G_er).astype(np.int8) print(A_ba.shape, A_ba.dtype)nx.to_numpy_array(G)按节点标签顺序返回 0/1 矩阵。N1000 时矩阵约 8MB作为 numpy 数组直接传进传播函数没有问题。.astype(np.int8)把 8 字节的 float64 压成 1 字节整数内存缩到八分之一。若节点标签不是从 0 开始的连续整数比如从 CSV 里读了“1、3、7”这种稀疏编号必须先重映射成 0…n-1否则矩阵的行列索引和节点对不上SIR 跑出来的邻居关系全是错的。这个坑在 5.4 和 6.3 还会再遇见。提示to_numpy_array在 networkx 2.x 之后的签名是to_numpy_array(G, nodelistNone, weightNone)默认按节点加入顺序排列。若要保证多次构建结果一致显式传入nodelistsorted(G.nodes())是最稳妥的。3. SIR 传播模块状态更新规则、同步更新与蒙特卡洛循环3.1 SIR 状态机与参数含义S→I→R 怎么流转SIR 模型把每个节点归入三种状态S易感、I感染、R恢复不再参与传播。每个离散时间步做两件事每个感染节点以概率 beta 尝试传染它的每个易感邻居同时每个感染节点以概率 gamma 转为恢复态。这个模型用两个参数就够beta 控制传染强度gamma 控制病程长度。当 beta 远大于 gamma 时疾病容易爆发当 beta 小到一定程度传播会随机死掉。工程上常用基础再生数 R0 beta / gamma 做第一把尺子但它只在均匀混合人群假设下严格成立。网络结构非均匀时R0 大于 1 并不保证爆发——BA 网络在 R0 略大于 1 时可能爆发ER 网络却可能静默。第四章会展开讲。同时要注意“同一时间步内”的先后顺序先传染再恢复还是先恢复再传染结果有细微差别。常见做法统一成“先传染后恢复”每个时间步结束再统计三种状态的数量。这样设定简单、可复现也符合大多数离散时间传播模型的习惯。3.2 最小可运行的单次传播函数以及同步更新设定下面这个函数是整套方案的核心单次传播直接靠它完成def sir_once(adj, beta, gamma, initial_inf, max_steps200, rngNone): n adj.shape[0] state np.zeros(n, dtypenp.int8) # 0: S, 1: I, 2: R state[initial_inf] 1 rng rng or np.random.default_rng() I_hist [int(np.count_nonzero(state 1))] R_hist [int(np.count_nonzero(state 2))] for _ in range(max_steps): infected np.flatnonzero(state 1) if infected.size 0: break # 1) 传染本轮在册的感染者去感染易感邻居 for i in infected: for j in np.flatnonzero(adj[i]): if state[j] 0 and rng.random() beta: state[j] 1 # 新感染者下一轮才具备传染性 # 2) 恢复本轮在册的感染者按概率 gamma 转 R for i in infected: if rng.random() gamma: state[i] 2 I_hist.append(int(np.count_nonzero(state 1))) R_hist.append(int(np.count_nonzero(state 2))) return np.array(I_hist), np.array(R_hist)逻辑说明infected在每轮开头取一次快照本轮新感染的节点不会立刻去传染别人而是等下一轮循环开始后再进入传播队列。这是标准的“同步更新”语义。同步更新的好处是结果与节点遍历顺序无关无论从哪个节点先开始遍历概率分布都一致代价是同样参数下传播速度比连续时间模型略慢。如果你要严格模拟连续时间过程就需要把感染和恢复事件按指数间隔排进事件队列逐事件推进这个方案里不展开。参数说明adj是上一章的 0/1 邻接矩阵beta是每对 S-I 邻居在每个时间步的传染概率gamma是感染者每时间步的恢复概率两者取值都在 0 到 1initial_inf是初始感染节点的下标列表max_steps是最大迭代轮数。返回的I_hist和R_hist第 0 个元素是 t0 时刻的初始状态。直观理解一个感染者平均存在 1/gamma 个时间步的病程每个时间步内向每个易感邻居发起一次 beta 概率的传染尝试。3.3 蒙特卡洛封装单次运行说明不了问题要跑上百次单次运行的随机性很大。同样的 beta、gamma 和网络跑三次可能得到三条完全不同的曲线小规模网络上尤其明显。我一般对同一组参数跑 100 到 200 次把所有循环的 I 和 R 曲线按时间步对齐后取平均。这里有一个高频踩坑点传播可能在第 80 步结束也可能拖到第 150 步直接把曲线叠成 numpy 数组会因为长度不一致而报错。做法是用最后一个值填充到最大长度——恢复数结束后本来就不变感染数结束后恒为 0填充不改变统计含义。def align_curve(curve, length): if len(curve) length: return curve[:length] return np.pad(curve, (0, length - len(curve)), modeedge) def sir_monte_carlo(adj, beta, gamma, initial_inf, runs100, max_steps200, seed42): rng np.random.default_rng(seed) all_I, all_R [], [] for _ in range(runs): I_hist, R_hist sir_once(adj, beta, gamma, initial_inf, max_stepsmax_steps, rngrng) all_I.append(I_hist) all_R.append(R_hist) length max(len(c) for c in all_I) I_mat np.vstack([align_curve(c, length) for c in all_I]) R_mat np.vstack([align_curve(c, length) for c in all_R]) return I_mat.mean(axis0), R_mat.mean(axis0)逻辑说明np.vstack把所有曲线堆成 runs 行 × length 列的矩阵axis0求平均得到平均感染人数曲线和平均恢复人数曲线。想要离散程度的话再把I_mat.std(axis0)一起返回。传播结束后感染人数为 0因此用modeedge填充感染曲线不会引入虚假数值恢复曲线则自然保持最终规模。种子在这里的作用很关键同一个 seed 传入default_rng整个蒙特卡洛过程的随机数序列就完全确定固定 seed 后别人复跑你的脚本能得到一模一样的曲线修改 seed 则可以检验结论对随机性的稳健性。种子一旦固定后续所有实验都基于同一份随机序列不同网络的对比不会因为某一次运气而偏移。4. 参数对照实验感染率、恢复率与初始感染节点怎么设4.1 基准参数表把实验条件固定下来做对比实验前先把一张参数表写清楚不然跑完三天再看数据可能已经想不起来当初的网络规模和初始感染是怎么设的。我常用的基准场景如下参数值说明节点数 N1000适中单次传播毫秒级BA 连边数 m3平均度约 6ER 连边概率 p2m/(N-1)期望平均度 6与 BA 对齐传染概率 beta0.2基准值实验中上下扫描恢复概率 gamma0.1平均病程 10 个时间步初始感染节点数10随机抽取固定 seed蒙特卡洛次数200保证统计稳定性最大时间步300确保传播结束或接近结束这张表里最值得解释的是 beta 和 gamma 的搭配。gamma0.1 意味着每个感染者平均在 10 步内恢复beta0.2 表示每个易感邻居每步有 20% 概率被传染。如果不把两者放在一起看很难判断 0.2 算强还是弱有了 gamma 做参照就能估算单个感染者链条的延续时间尺度。4.2 扫描 beta对比最终感染规模固定 gamma 和初始感染者后把 beta 从 0.02 扫到 0.5统计每个 beta 下的最终恢复比例 rho R[-1]/N。rho 就是“疫情最终波及多大比例的网络”是传播动力学中最常用的指标之一。def epidemic_size(adj, beta, gamma, initial_inf, runs200, seed42): _, R_hist sir_monte_carlo(adj, beta, gamma, initial_inf, runsruns, seedseed) return R_hist[-1] / adj.shape[0] betas np.linspace(0.02, 0.5, 15) rho_ba [epidemic_size(A_ba, b, 0.1, initial_inf, seed42) for b in betas] rho_er [epidemic_size(A_er, b, 0.1, initial_inf, seed42) for b in betas] for b, rb, re in zip(betas, rho_ba, rho_er): print(fbeta{b:.3f} BA{rb:.3f} ER{re:.3f})逻辑说明sir_monte_carlo内部每次从 seed42 重新生成随机数序列所以任何 beta 取值下的 200 次重复都完全可复现扫描过程中唯一的变量就是 beta 本身。R_hist的最后一个元素是传播稳定后处于恢复态的节点数因为恢复态不可逆它等于“曾经被感染过的人数”除以 N 就是最终感染规模。典型的实验结果beta 很小时比如 0.05ER 网络几乎不爆发rho 只有百分之几BA 网络在同样的 beta 下 rho 已经明显抬升。等到 beta 超过 0.3差距又会缩小因为强传染参数下两个网络都会被扫一遍结构差异被淹没。如果你跑出来的曲线看不出这个趋势第一件事检查网络是否对齐第二件事确认两张图的initial_inf用的是不是同一个列表。需要注意的是 beta 和 gamma 的另一种扫描方式。若把 gamma 从 0.1 改到 0.5平均病程从 10 步缩到 2 步同一个感染者在网络里停留时间变短传给邻居的机会骤减。gamma 变化不单纯是“恢复快慢”它还直接改变了基础再生数 beta/gamma。所以入门阶段最干净的实验设计是固定 gamma 只扫 beta如果一定要扫 gamma把横坐标换成 beta/gamma 后再对比两个网络两条曲线的结构差异才会从噪声中凸显出来。4.3 BA 网络传播阈值趋零的机制hub 节点和度方差为什么结构差异会带来阈值差异先给结论再解释。常用近似是异质平均场它把传播阈值写成 beta_c ⟨k⟩ / ⟨k²⟩。ER 网络度分布接近泊松⟨k²⟩ ≈ ⟨k⟩² ⟨k⟩代进去得到 beta_c ≈ 1/⟨k⟩在 m3 时约 0.167。BA 网络的度分布是幂律⟨k²⟩ 被大度节点显著抬高有限网络里算出来通常只有 0.02 到 0.05理论极限下趋近于 0。通俗地讲ER 网络里大家度都差不多一个感染者很难同时面对大量易感邻居BA 网络里只要一个 hub 被传染它会同时面对几十上百个易感邻居发起传染尝试即使单次传染概率很低总有一次会成功。hub 把传播链“接通”的能力让 BA 网络在 beta 极低时仍然可能形成大规模流行。这同时也是 5.5 里“初始感染者选 hub 还是随机节点会改变结果”的理论根源。要补充一句这个阈值只是定性参考。有限尺寸的随机网络在阈值附近震荡非常剧烈1000 节点的单次实验结果在阈值附近跳跃性很大。如果你的目标是计算精确阈值需要把网络尺寸拉到 5000 以上并做大量重复目标是对比两个网络结构差异导致的行为差异1000 节点完全够用。我做对比实验时还会额外看一个指标平均感染峰值出现的时间步。ER 网络通常要更晚才到峰值因为传播需要逐步扩散到整个网络BA 网络因为有 hub 作为捷径峰值来得更早。这个时序差异在论文图表里非常直观适合和最终规模一起展示。5. 避坑清单五个让 BA/ER 对比翻车的常见问题5.1 随机种子不固定结果玄学化现象同一份代码今天跑出一个流行曲线明天跑出另一个谁都不敢采信。原因np.random.random()在没有显式种子时从系统熵源取随机数每次运行序列都不一样。SIR 过程又极端依赖随机序列小网络里完全可能第一次爆发、第二次绝灭单次结果之间的方差大到足以淹没结构差异。解决不要依赖全局随机状态。在sir_once和sir_monte_carlo里显式创建np.random.default_rng(seed)把 rng 对象作为参数层层传递。记录实验时把“随机种子 初始感染节点 网络参数”三件套写进注释或文件名三个月后回来看还能复现。额外注意一点如果循环结构变动导致随机数调用次数变化即使 seed 相同后续序列也会整体偏移。因此最好统一让 rng 对象只被传播函数消费上层代码不要在中间插入其他随机数调用。5.2 平均度或网络规模不一致比较失去意义现象BA 和 ER 的参数随便设跑出来 ER 传播得又快又广于是得出“ER 更容易传播”的错误结论。原因ER 的 p 给大了边数暴涨传播路径成倍增加或者两个网络的节点数 n 不统一大网络本身更容易容留传播链。对比实验的前提是“唯一变量是网络结构”边数和规模必须一致。解决生成两张图后立刻打印number_of_nodes、number_of_edges和平均度用p 2m/(n-1)对齐。我习惯在实验脚本开头加断言检查两个网络边数差是否在 5% 以内防止之后改动参数时漏掉同步修改。如果边数差异很大先回去检查 p 的推导而不是急着分析传播曲线。5.3 恢复概率 gamma 与病程的关系没算清现象把 gamma 从 0.1 改成 0.5然后和 gamma0.1 的曲线对比得出“传染率变高所以传播更快”的错误结论。原因gamma0.5 意味着平均病程只有 2 步感染者在网络里停留时间短传给邻居的机会大幅减少。gamma 同时影响“病程”和“基础再生数 beta/gamma”不是一个单纯的速度旋钮。解决固定 gamma 只扫 beta是入门阶段最干净的实验设计。要扫 gamma 时把横坐标改为 beta/gamma基础再生数两个网络各自的结果才能放在同一尺度下比较。但在低 beta 区域这种重标定也会掩盖阈值形状差异所以我建议至少加一组“固定 gamma 扫 beta”的对照确保结构差异没有被参数重标定抹掉。5.4 邻接矩阵全量存储网络一大就卡死现象把 N 从 1000 改成 10000nx.to_numpy_array生成 float64 的 10000×10000 矩阵瞬间占掉约 800MB 内存传播循环慢到无法忍受。原因稠密邻接矩阵的存储是 O(n²)10000²×8 字节就是 800MB加上sir_once里两层 Python 循环遍历邻居每步都重扫矩阵行慢是必然的。解决先用.astype(np.int8)把 0/1 矩阵压到八分之一内存。对 n5000 的网络改用scipy.sparse.csr_matrixsir_once里用indices和indptr按行取邻居避免np.flatnonzero(adj[i])扫描整行。判断标准1000 节点以下用 int8 稠密矩阵足够要到 10 万节点就必须把传播逻辑迁到稀疏表示这一步省不掉。5.5 初始感染节点选 hub 还是随机节点结果差异显著现象把初始感染者从随机 10 个改成“度最大的 10 个”BA 网络的最终规模从 30% 跳到 80%而 ER 网络几乎不变。原因BA 网络的 hub 是传播的超级中转站。初始感染直接从 hub 开始等于跳过了流行形成最困难的随机阶段ER 网络没有这样的 hub初始位置的影响小得多。这个效应不是 bug而是结构特性但如果不控制它会成为实验里一个巨大的混淆变量。解决对比实验里必须统一初始感染策略。常见做法是随机抽 K 个节点用固定 seed 生成那组下标两个网络使用完全相同的初始感染者列表如果故意研究初始位置的影响把“随机 vs 最大度”拆成单独一组实验。写报告时务必在方法部分写清“初始感染为随机抽样且种子固定”否则别人无法复现结果。6. 验证与进阶从理论阈值对照到自有数据迁移6.1 把仿真结果和理论阈值放在一张图里仿真跑完先别急着下结论。把两个网络的实际度序列拿出来算 ⟨k⟩ 和 ⟨k²⟩得到各自的异质平均场预测阈值for name, G in [(BA, G_ba), (ER, G_er)]: k np.array([d for _, d in G.degree()]) kbar k.mean() k2bar (k**2).mean() print(f{name}: k{kbar:.2f} k^2{k2bar:.1f} threshold{kbar/k2bar:.4f})把算出的阈值竖线画在 rho-beta 图上应该看到BA 的竖线落在很靠左的位置ER 的竖线明显偏右仿真曲线抬升的顺序也应该与之一致。如果出现 BA 和 ER 的抬升次序反过来的情况优先回头查代码而不是怀疑理论——这通常意味着平均度没对齐或者初始感染节点选取不一致。6.2 把传播过程存成动画与可读图光看数字不够直观。把 I 曲线和 R 曲线画成折线图保存横坐标是时间步时间步超过 300 时刻度会挤成一团用plt.xticks(np.arange(0, len(I_hist), 30))手动指定步长。想观察网络中哪些区域先被感染可以用nx.spring_layout布局 1000 节点的网络配合FuncAnimation把每步的感染节点标红存成 GIF节点到 3000 以上就放弃布局动画改拍感染人数曲线动画更实际。import matplotlib.pyplot as plt plt.figure(figsize(8, 5)) plt.plot(betas, rho_ba, markero, labelBA) plt.plot(betas, rho_er, markers, labelER) plt.xlabel(beta) plt.ylabel(final epidemic size) plt.legend() plt.tight_layout() plt.savefig(rho_beta.png, dpi150)这样保存图片并调整标注密度比在终端里看一串数字直观得多也方便直接贴进实验报告。dpi150足够印刷需求文件也不会太大。6.3 从 BA/ER 迁移到自有网络数据的三个改动真实场景里你手里往往是一张自己的网络而不是现成的 BA/ER 图。迁移时要改三处第一从 CSV 或 txt 读边列表节点编号先重映射到 0 到 n-1再用nx.from_edgelist建图第二确认网络是否有向——有向图要把“对所有邻居”改成“对所有出边邻居”加权网络则把边的权重乘到 beta 上第三检查自环和多重边自环对传播没有意义多重边会让邻接矩阵不再是 0/1建图时必须先去重否则度分布统计会被顶歪。我的个人习惯是每次实验跑完至少留三样东西——随机种子、初始感染节点列表、网络参数表。最初我对比 BA 和 ER 时初始感染节点随手随机抽了两次结果 ER 一度比 BA 传播得还猛排查到凌晨才发现是初始节点位置不同造成的。从那以后种子和参数表成了我所有仿真项目的固定开头。这个教训让我明白传播动力学实验里 90% 的“惊喜”其实都是变量没控制住。希望这篇笔记帮你绕开这些踩过的坑把时间真正花在看清楚网络结构如何塑造传播动力学上。本文还有配套的精品资源点击获取