向前追踪与向后追踪——传播树中的朋友悖论

# 在度异质网络上做 SIR 模拟并记录“谁传染了谁”（传播树）。
# 比较：随机挑一个感染者，(a) 它自己的继发病例数（向前追踪能发现的人数），
#       (b) 它的“传染源”的继发病例数（向后追踪找到的源头还能牵出多少人）。
# 解读：整个暴发中“每个病例平均只传给约 1 人”（除种子外每个病例恰有一个传染源），
#       但随机病例的传染源平均传给了多得多的人——被感染者更可能来自“多产”的传染源，
#       这正是 Kojaku 等 2021 指出向后追踪更有效的原因。
import numpy as np
import networkx as nx

rng = np.random.default_rng(5)
G = nx.barabasi_albert_graph(5000, 2, seed=5)
beta, gamma = 0.15, 0.3

def sir_tree(G, seed, rng):
    parent = {seed: None}; state = {seed: 1}; frontier = [seed]
    while frontier:
        nxt = []
        for u in frontier:
            for v in G.neighbors(u):
                if v not in state and rng.random() < beta:
                    state[v] = 1; parent[v] = u; nxt.append(v)
            if rng.random() < gamma:
                state[u] = 2
            else:
                nxt.append(u)
        frontier = [v for v in nxt if state[v] == 1]
    return parent

own, src_off = [], []
for _ in range(20):
    parent = sir_tree(G, int(rng.integers(5000)), rng)
    if len(parent) < 200:
        continue
    off = {}
    for v, p in parent.items():
        if p is not None:
            off[p] = off.get(p, 0) + 1
    for v, p in parent.items():
        own.append(off.get(v, 0))
        if p is not None:
            src_off.append(off[p])
print("随机感染者的平均继发病例数（≈R）= %.2f" % np.mean(own))
print("其传染源的平均继发病例数       = %.2f" % np.mean(src_off))
print("继发病例数的方差/均值 = %.2f（>1 表示过度离散）" % (np.var(own) / np.mean(own)))

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 4.2))
bins = np.arange(0, max(max(own), max(src_off)) + 2) - 0.5
ax.hist(own, bins=bins, density=True, histtype="step", lw=2, label="random case (forward)")
ax.hist(src_off, bins=bins, density=True, histtype="step", lw=2, ls="--", label="its infector (backward)")
ax.axvline(np.mean(own), c="tab:blue", lw=0.8); ax.axvline(np.mean(src_off), c="tab:orange", lw=0.8)
ax.set_yscale("log"); ax.set_xlabel("number of secondary cases"); ax.set_ylabel("probability"); ax.legend()
ax.set_title("Infectors have more offspring (means %.1f vs %.1f)" % (np.mean(own), np.mean(src_off)))
plt.tight_layout(); plt.show()