Example05 / 14

三种阈值近似的比较——异质平均场、淬火平均场与非回溯矩阵

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

同一张网络可以给出三个不同的 SIR 阈值估计： HMF：T_c = /(<k^2>-)；QMF：T_c ≈ 1/Λ_max(A)；消息传递：T_c = 1/ρ(B)，B 为非回溯矩阵。 ρ(B) 用 Ihara–Bass 公式化成 2N×2N 矩阵 [[A, I-D], [I, 0]] 的最大特征值来算，避免构造 2m×2m 的 B。再用“键渗流”精确地模拟 SIR 最终规模（每条边以概率 T 保留，取种子所在分支），看哪一个估计最接近实际转变点。解读：三者的大小关系一般为 QMF ≤ 非回溯，且 1/Λ_max 常偏低（Λ_max ≥ sqrt(k_max)，被单个枢纽拉高）；有限网络中转变是平滑的，最大分支在非回溯估计附近开始明显离开 0。

Output输出图

Output文字输出

N = 3000, <k> = 4.90, k_max = 240
HMF  T_c = 0.0505
QMF  T_c = 1/Λ_max = 0.0481
非回溯 T_c = 1/ρ(B) = 0.0612
T = 0.03  最大分支（大暴发规模）= 0.0048
T = 0.05  最大分支（大暴发规模）= 0.0134
T = 0.07  最大分支（大暴发规模）= 0.0353
T = 0.09  最大分支（大暴发规模）= 0.0617
T = 0.12  最大分支（大暴发规模）= 0.1297
T = 0.15  最大分支（大暴发规模）= 0.2024
T = 0.20  最大分支（大暴发规模）= 0.3162

Code代码

# 同一张网络可以给出三个不同的 SIR 阈值估计：
#   HMF：T_c = <k>/(<k^2>-<k>)；QMF：T_c ≈ 1/Λ_max(A)；消息传递：T_c = 1/ρ(B)，B 为非回溯矩阵。
# ρ(B) 用 Ihara–Bass 公式化成 2N×2N 矩阵 [[A, I-D], [I, 0]] 的最大特征值来算，避免构造 2m×2m 的 B。
# 再用“键渗流”精确地模拟 SIR 最终规模（每条边以概率 T 保留，取种子所在分支），看哪一个估计最接近实际转变点。
# 解读：三者的大小关系一般为 QMF ≤ 非回溯，且 1/Λ_max 常偏低（Λ_max ≥ sqrt(k_max)，被单个枢纽拉高）；
#       有限网络中转变是平滑的，最大分支在非回溯估计附近开始明显离开 0。
import numpy as np
import networkx as nx
import scipy.sparse as sp
import scipy.sparse.linalg as spla

rng = np.random.default_rng(2)
deg_seq = np.clip(np.round(rng.pareto(1.6, 3000) * 2 + 2), 2, 300).astype(int)
if deg_seq.sum() % 2:
    deg_seq[0] += 1
G = nx.Graph(nx.configuration_model(deg_seq, seed=2))
G.remove_edges_from(nx.selfloop_edges(G))
G = G.subgraph(max(nx.connected_components(G), key=len)).copy()
G = nx.convert_node_labels_to_integers(G)
N = G.number_of_nodes()

k = np.array([d for _, d in G.degree()], dtype=float)
Tc_hmf = k.mean() / ((k**2).mean() - k.mean())
A = nx.to_scipy_sparse_array(G, format="csr", dtype=float)
lam_max = spla.eigsh(A, k=1, which="LA")[0][0]
I = sp.identity(N, format="csr")
D = sp.diags(k)
M = sp.bmat([[A, I - D], [I, None]], format="csr")
rho_B = np.max(np.abs(spla.eigs(M, k=1, which="LM")[0]))
print("N = %d, <k> = %.2f, k_max = %d" % (N, k.mean(), k.max()))
print("HMF  T_c = %.4f" % Tc_hmf)
print("QMF  T_c = 1/Λ_max = %.4f" % (1 / lam_max))
print("非回溯 T_c = 1/ρ(B) = %.4f" % (1 / rho_B))

def giant_after_percolation(G, T, rng, runs=30):
    """键渗流：每条边以概率 T 保留，返回最大连通分支的平均相对大小（即大暴发的最终规模）。"""
    edges = np.array(G.edges())
    out = []
    for _ in range(runs):
        H = nx.Graph(); H.add_nodes_from(G); H.add_edges_from(edges[rng.random(len(edges)) < T])
        out.append(max(len(c) for c in nx.connected_components(H)) / N)
    return np.mean(out)

Ts, giants = [0.03, 0.05, 0.07, 0.09, 0.12, 0.15, 0.20], []
for T in Ts:
    giants.append(giant_after_percolation(G, T, rng))
    print("T = %.2f  最大分支（大暴发规模）= %.4f" % (T, giants[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.plot(Ts, giants, "ko-", label="bond percolation (simulated)")
for x, name, ls in [(Tc_hmf, "HMF", ":"), (1 / lam_max, "QMF 1/Lambda_max", "--"), (1 / rho_B, "non-backtracking 1/rho(B)", "-.")]:
    ax.axvline(x, ls=ls, c="tab:red" if ls == "-." else "0.4", label="%s = %.3f" % (name, x))
ax.set_xlabel("transmissibility T"); ax.set_ylabel("largest outbreak size / N"); ax.legend(fontsize=8)
ax.set_title("Outbreak size departs from 0 near the NB threshold")
plt.tight_layout(); plt.show()

Run在浏览器中运行

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

三种阈值近似的比较——异质平均场、淬火平均场与非回溯矩阵

# 同一张网络可以给出三个不同的 SIR 阈值估计：
#   HMF：T_c = <k>/(<k^2>-<k>)；QMF：T_c ≈ 1/Λ_max(A)；消息传递：T_c = 1/ρ(B)，B 为非回溯矩阵。
# ρ(B) 用 Ihara–Bass 公式化成 2N×2N 矩阵 [[A, I-D], [I, 0]] 的最大特征值来算，避免构造 2m×2m 的 B。
# 再用“键渗流”精确地模拟 SIR 最终规模（每条边以概率 T 保留，取种子所在分支），看哪一个估计最接近实际转变点。
# 解读：三者的大小关系一般为 QMF ≤ 非回溯，且 1/Λ_max 常偏低（Λ_max ≥ sqrt(k_max)，被单个枢纽拉高）；
#       有限网络中转变是平滑的，最大分支在非回溯估计附近开始明显离开 0。
import numpy as np
import networkx as nx
import scipy.sparse as sp
import scipy.sparse.linalg as spla

rng = np.random.default_rng(2)
deg_seq = np.clip(np.round(rng.pareto(1.6, 3000) * 2 + 2), 2, 300).astype(int)
if deg_seq.sum() % 2:
    deg_seq[0] += 1
G = nx.Graph(nx.configuration_model(deg_seq, seed=2))
G.remove_edges_from(nx.selfloop_edges(G))
G = G.subgraph(max(nx.connected_components(G), key=len)).copy()
G = nx.convert_node_labels_to_integers(G)
N = G.number_of_nodes()

k = np.array([d for _, d in G.degree()], dtype=float)
Tc_hmf = k.mean() / ((k**2).mean() - k.mean())
A = nx.to_scipy_sparse_array(G, format="csr", dtype=float)
lam_max = spla.eigsh(A, k=1, which="LA")[0][0]
I = sp.identity(N, format="csr")
D = sp.diags(k)
M = sp.bmat([[A, I - D], [I, None]], format="csr")
rho_B = np.max(np.abs(spla.eigs(M, k=1, which="LM")[0]))
print("N = %d, <k> = %.2f, k_max = %d" % (N, k.mean(), k.max()))
print("HMF  T_c = %.4f" % Tc_hmf)
print("QMF  T_c = 1/Λ_max = %.4f" % (1 / lam_max))
print("非回溯 T_c = 1/ρ(B) = %.4f" % (1 / rho_B))

def giant_after_percolation(G, T, rng, runs=30):
    """键渗流：每条边以概率 T 保留，返回最大连通分支的平均相对大小（即大暴发的最终规模）。"""
    edges = np.array(G.edges())
    out = []
    for _ in range(runs):
        H = nx.Graph(); H.add_nodes_from(G); H.add_edges_from(edges[rng.random(len(edges)) < T])
        out.append(max(len(c) for c in nx.connected_components(H)) / N)
    return np.mean(out)

Ts, giants = [0.03, 0.05, 0.07, 0.09, 0.12, 0.15, 0.20], []
for T in Ts:
    giants.append(giant_after_percolation(G, T, rng))
    print("T = %.2f  最大分支（大暴发规模）= %.4f" % (T, giants[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.plot(Ts, giants, "ko-", label="bond percolation (simulated)")
for x, name, ls in [(Tc_hmf, "HMF", ":"), (1 / lam_max, "QMF 1/Lambda_max", "--"), (1 / rho_B, "non-backtracking 1/rho(B)", "-.")]:
    ax.axvline(x, ls=ls, c="tab:red" if ls == "-." else "0.4", label="%s = %.3f" % (name, x))
ax.set_xlabel("transmissibility T"); ax.set_ylabel("largest outbreak size / N"); ax.legend(fontsize=8)
ax.set_title("Outbreak size departs from 0 near the NB threshold")
plt.tight_layout(); plt.show()