Example02 / 14

用度分布估计 SIR 流行病阈值

所属专题: 网络上的传播与免疫
下载: Notebook · .py

异质平均场给出 SIR 的阈值：传播性 T 需大于 T_c = / (<k^2> - )。在上一示例的离散时间模型中，一个感染者在整个感染期内把病传给某个邻居的概率（传播性）为 T = 1 - gamma(1-beta) / (1 - (1-beta)(1-gamma)) 代码计算 ER 与 BA 网络的 T_c，并与当前参数的 T 比较。解读：BA 网络的 <k^2> 更大，T_c 更低；同一 T 可能高于 BA 的阈值而低于 ER 的阈值。

Output输出图

Output文字输出

当前参数的传播性 T = 0.149
ER: T_c = 0.167  -> 低于阈值
BA: T_c = 0.082  -> 高于阈值，可能暴发

Code代码

# 异质平均场给出 SIR 的阈值：传播性 T 需大于 T_c = <k> / (<k^2> - <k>)。
# 在上一示例的离散时间模型中，一个感染者在整个感染期内把病传给某个邻居的概率（传播性）为
#   T = 1 - gamma*(1-beta) / (1 - (1-beta)*(1-gamma))
# 代码计算 ER 与 BA 网络的 T_c，并与当前参数的 T 比较。
# 解读：BA 网络的 <k^2> 更大，T_c 更低；同一 T 可能高于 BA 的阈值而低于 ER 的阈值。
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

N, m = 1000, 3
G_ba = nx.barabasi_albert_graph(N, m, seed=1)
G_er = nx.gnm_random_graph(N, G_ba.number_of_edges(), seed=1)
beta, gamma = 0.05, 0.3
T = 1 - gamma * (1 - beta) / (1 - (1 - beta) * (1 - gamma))
print("当前参数的传播性 T = %.3f" % T)

def threshold(G):
    k = np.array([d for _, d in G.degree()], dtype=float)
    return k.mean() / ((k**2).mean() - k.mean())

for name, G in [("ER", G_er), ("BA", G_ba)]:
    Tc = threshold(G)
    print("%s: T_c = %.3f  -> %s" % (name, Tc, "高于阈值，可能暴发" if T > Tc else "低于阈值"))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
for name, G, mk in [("ER", G_er, "o"), ("BA", G_ba, "s")]:
    k, cnt = np.unique([d for _, d in G.degree()], return_counts=True)
    ax1.loglog(k, cnt / N, mk, ms=4, label=name)
ax1.set_xlabel("degree k"); ax1.set_ylabel("P(k)"); ax1.legend(); ax1.set_title("Degree distributions")
ax2.bar(["ER", "BA"], [threshold(G_er), threshold(G_ba)], color=["0.6", "tab:blue"], width=0.5)
ax2.axhline(T, ls="--", c="tab:red", label="T = %.3f (current)" % T)
ax2.set_ylabel("epidemic threshold T_c"); ax2.legend()
ax2.set_title("Heterogeneity lowers T_c: BA outbreak, ER not")
plt.tight_layout(); plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

用度分布估计 SIR 流行病阈值

# 异质平均场给出 SIR 的阈值：传播性 T 需大于 T_c = <k> / (<k^2> - <k>)。
# 在上一示例的离散时间模型中，一个感染者在整个感染期内把病传给某个邻居的概率（传播性）为
#   T = 1 - gamma*(1-beta) / (1 - (1-beta)*(1-gamma))
# 代码计算 ER 与 BA 网络的 T_c，并与当前参数的 T 比较。
# 解读：BA 网络的 <k^2> 更大，T_c 更低；同一 T 可能高于 BA 的阈值而低于 ER 的阈值。
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

N, m = 1000, 3
G_ba = nx.barabasi_albert_graph(N, m, seed=1)
G_er = nx.gnm_random_graph(N, G_ba.number_of_edges(), seed=1)
beta, gamma = 0.05, 0.3
T = 1 - gamma * (1 - beta) / (1 - (1 - beta) * (1 - gamma))
print("当前参数的传播性 T = %.3f" % T)

def threshold(G):
    k = np.array([d for _, d in G.degree()], dtype=float)
    return k.mean() / ((k**2).mean() - k.mean())

for name, G in [("ER", G_er), ("BA", G_ba)]:
    Tc = threshold(G)
    print("%s: T_c = %.3f  -> %s" % (name, Tc, "高于阈值，可能暴发" if T > Tc else "低于阈值"))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
for name, G, mk in [("ER", G_er, "o"), ("BA", G_ba, "s")]:
    k, cnt = np.unique([d for _, d in G.degree()], return_counts=True)
    ax1.loglog(k, cnt / N, mk, ms=4, label=name)
ax1.set_xlabel("degree k"); ax1.set_ylabel("P(k)"); ax1.legend(); ax1.set_title("Degree distributions")
ax2.bar(["ER", "BA"], [threshold(G_er), threshold(G_ba)], color=["0.6", "tab:blue"], width=0.5)
ax2.axhline(T, ls="--", c="tab:red", label="T = %.3f (current)" % T)
ax2.set_ylabel("epidemic threshold T_c"); ax2.legend()
ax2.set_title("Heterogeneity lowers T_c: BA outbreak, ER not")
plt.tight_layout(); plt.show()