数学推导与算法

本页补全前面各页用到的结果：A. SIS 阈值 $\lambda_c=\langle k\rangle/\langle k^2\rangle$ 的推导；B. SIR 最终规模与键渗流的对应；C. 精确的随机模拟算法；D. 活动驱动时间网络上的阈值；E. 淬火平均场与非回溯矩阵阈值；F. 目标免疫所需的接种比例；G. 阈值模型的全局级联条件。记号：$P(k)$ 为度分布，$\langle k^n\rangle=\sum_k k^nP(k)$，$\beta$ 为每条 S–I 边上的传播速率，$\mu$ 为恢复速率，$\lambda=\beta/\mu$。

A. SIS 的异质平均场阈值

假设（近似在这里进入）：（1）网络是无度相关的配置模型，$N\to\infty$；（2）同一度类的节点在统计上等价，只用 $i_k$ 描述（“退火”近似：相当于每一时刻重新随机连边）；（3）忽略邻居状态之间的相关性。

第 1 步：$\Theta$ 的含义。 从任一节点出发沿一条随机边走，到达度为 $k'$ 的节点的概率与它拥有的“边端”数成正比，为 $k'P(k')/\langle k\rangle$。因此该邻居处于感染态的概率为 $$\Theta=\frac{1}{\langle k\rangle}\sum_{k'}k'P(k')\,i_{k'} .$$

第 2 步：稳态。 令 HMF 方程 $\dfrac{di_k}{dt}=\beta k(1-i_k)\Theta-\mu i_k$ 等于零，两边除以 $\mu$： $$i_k=\frac{\lambda k\Theta}{1+\lambda k\Theta}.$$

第 3 步：自洽方程。 把 $i_k$ 代回 $\Theta$ 的定义： $$\Theta=F(\Theta)\equiv\frac{1}{\langle k\rangle}\sum_k kP(k)\,\frac{\lambda k\Theta}{1+\lambda k\Theta}.$$ $\Theta=0$ 总是解（无病态）。$F$ 满足 $F(0)=0$、$F(1)<1$，且 $F''(\Theta)<0$（每一项都是 $\Theta$ 的凹函数）。因此非零解存在，当且仅当 $F$ 在原点的斜率大于 1： $$F'(0)=\frac{\lambda}{\langle k\rangle}\sum_k k^2P(k)=\lambda\frac{\langle k^2\rangle}{\langle k\rangle}>1 \quad\Longleftrightarrow\quad \lambda>\lambda_c=\frac{\langle k\rangle}{\langle k^2\rangle}.$$

第 4 步：为什么 $\gamma\le3$ 时阈值消失。 取 $P(k)=Ck^{-\gamma}$，$k_{\min}\le k\le k_{\max}$，用连续近似： $$\langle k^2\rangle\approx C\int_{k_{\min}}^{k_{\max}}k^{2-\gamma}\,dk= \begin{cases} \dfrac{C\,(k_{\max}^{3-\gamma}-k_{\min}^{3-\gamma})}{3-\gamma}, & \gamma\ne3,\\[2mm] C\ln(k_{\max}/k_{\min}), & \gamma=3 . \end{cases}$$ $2<\gamma<3$ 时指数 $3-\gamma>0$，$\langle k^2\rangle$ 随 $k_{\max}$ 发散；$\gamma=3$ 时对数发散；$\gamma>3$ 时收敛到有限值。配置模型中“自然截断”$k_{\max}\sim k_{\min}N^{1/(\gamma-1)}$，代入得 $2<\gamma<3$ 时 $\lambda_c\sim N^{-(3-\gamma)/(\gamma-1)}$，$\gamma=3$ 时 $\lambda_c\sim1/\ln N$。阈值随 $N$ 增大而趋于零，但下降很慢——这就是“有限网络阈值小而非零”的定量含义。

第 5 步：早期增长（SI 的特征时间）。 对 SI（$\mu=0$）在 $i_k\ll1$ 时线性化，并注意一个被感染节点已经用掉一条边（感染它的那条），定义 $\Theta=\sum_k(k-1)P(k)i_k/\langle k\rangle$，则 $d\Theta/dt=\beta\,\dfrac{\langle k^2\rangle-\langle k\rangle}{\langle k\rangle}\,\Theta$，得到度异质性与流行病阈值表中的 $\tau=\langle k\rangle/[\beta(\langle k^2\rangle-\langle k\rangle)]$。SIS 与 SIR 的特征时间按同样方式加入恢复项得到。

B. SIR 最终规模与键渗流

第 1 步：传播性。 设感染者 $i$ 的感染期为 $\tau$，在此期间沿每条边以速率 $\beta$ 传播。邻居 $j$ 最终被 $i$ 传染（假设 $j$ 尚未被别人传染）的概率为 $T=1-e^{-\beta\tau}$。若感染期服从速率为 $\mu$ 的指数分布，对 $\tau$ 求平均得 $T=\int_0^\infty\mu e^{-\mu\tau}(1-e^{-\beta\tau})d\tau=\dfrac{\beta}{\beta+\mu}$。

第 2 步：映射。 若 $\tau$ 对所有人相同，每条边是否“传得过去”是独立的伯努利事件，概率为 $T$。于是“最终被感染的集合”恰好等于“以概率 $T$ 保留每条边后、包含初始感染者的连通分支”——这就是键渗流（bond percolation）。若 $\tau$ 因人而异，同一感染者的各条出边不再独立（它们共享同一个 $\tau$），映射只是部分成立（Kenah 与 Robins 2007，Physical Review E）：在树状近似下，用平均传播性 $T$ 算出的大爆发规模仍然正确，但发生大爆发的概率不再等于这一规模，一般会更小。

第 3 步：生成函数求解。 记 $G_0(x)=\sum_kP(k)x^k$，$G_1(x)=G_0'(x)/G_0'(1)$（沿边到达节点的“剩余度”生成函数）。令 $v$ 为“沿一条边出发、不能通向大爆发”的概率：这条边要么没传过去（概率 $1-T$），要么传过去了但对方的其余边都不通向大爆发（概率 $T\,G_1(v)$）： $$v=1-T+T\,G_1(v),\qquad R_\infty=1-G_0(v).$$ $R_\infty$ 是大爆发中最终被感染（移除）的比例，适用于初始感染者占比趋于零的极限。

第 4 步：阈值。 $v=1$ 总是解。与 A 中同理，存在 $v<1$ 的解当且仅当 $T\,G_1'(1)>1$，而 $G_1'(1)=\dfrac{\langle k^2\rangle-\langle k\rangle}{\langle k\rangle}$，故 $$T_c=\frac{\langle k\rangle}{\langle k^2\rangle-\langle k\rangle},$$ 这就是度异质性与流行病阈值中的 $\lambda_c^{\text{SIR}}$（那里的 $\lambda$ 应理解为传播性 $T$）。

检验：泊松网络 $G_0(x)=G_1(x)=e^{\langle k\rangle(x-1)}$。由 $R_\infty=1-G_0(v)$ 得 $G_1(v)=1-R_\infty$，代入 $v$ 的方程得 $v=1-TR_\infty$，于是 $$R_\infty=1-e^{-\langle k\rangle T R_\infty},$$ 正是均匀混合下的经典最终规模方程，其中 $R_0=\langle k\rangle T$。

与消息传递、EBCM 的关系：把第 3 步中“对所有边相同的 $v$”换成每条有向边各自的 $u_{i\leftarrow j}$，就得到配对近似、边基模型与消息传递中的消息传递方程；EBCM 中 $\theta(t\to\infty)$ 满足的方程在指数感染期下与第 3 步等价。

C. 精确随机模拟：Gillespie 算法

微分方程描述的是期望值；有限网络上的随机涨落（例如疫情早期灭绝）需要随机模拟。对指数分布的传播和恢复时间，Gillespie 算法给出与连续时间马尔可夫过程分布完全一致的轨迹。

算法  网络 SIR 的 Gillespie 模拟
输入：网络 G、β、μ、初始感染集 I0
维护：状态 x_i∈{S,I,R}；对每个 i∈I 记 s_i = i 的易感邻居数；N_SI = Σ_{i∈I} s_i
1. t ← 0；按 I0 初始化 x、s、N_SI
2. while I 非空:
3.     Λ ← μ|I| + β·N_SI                     # 总事件速率
4.     抽 r1, r2 ~ U(0,1)；t ← t − ln(r1)/Λ    # 等待时间 ~ Exp(Λ)
5.     if r2·Λ < μ|I|:                          # 恢复事件
6.         从 I 中均匀抽 i；x_i ← R；N_SI ← N_SI − s_i；从 I 删除 i
7.     else:                                    # 感染事件
8.         以 s_i 为权重从 I 中抽 i；在 i 的易感邻居中均匀抽 j
9.         x_j ← I；s_j ← 0
10.        for j 的每个邻居 l:
11.            if x_l = I: s_l ← s_l − 1；N_SI ← N_SI − 1
12.            if x_l = S: s_j ← s_j + 1；N_SI ← N_SI + 1
13.        把 j 加入 I
14.    记录 (t, |S|, |I|, |R|)

复杂度：SIR 中每个节点至多被感染一次、恢复一次，事件总数 $\le 2N$；第 10–12 步对每个被感染节点只遍历一次邻居，合计 $O(m)$；第 8 步按权重抽样若用二叉求和树实现为 $O(\log N)$，总计 $O(m+N\log N)$。若改为 SIS（第 6 步令 $x_i\leftarrow S$ 并更新其邻居的 $s$），事件数不再有界，需设定最长模拟时间。离散时间近似（每步 $\Delta t$ 内以 $1-e^{-\beta\Delta t}$ 传播、$1-e^{-\mu\Delta t}$ 恢复）实现更简单，但有 $O(\Delta t)$ 的离散化误差。

事件驱动的等价算法。对 SIR 还有一种更简单的精确方法（Kiss、Miller 与 Simon 2017 书中称为快速 SIR 算法）：节点 $v$ 在时刻 $t$ 被感染时，抽取恢复时刻 $t+\mathrm{Exp}(\mu)$，并为每个邻居抽取传播时刻 $t+\mathrm{Exp}(\beta)$；只有早于恢复时刻的传播事件进入优先队列。按时间顺序处理队列，跳过已非易感的目标。代码示例“边基仓室模型与精确随机模拟的比较”使用的就是这一算法。把 $\mathrm{Exp}$ 换成其他分布，就能模拟非马尔可夫的感染期（例如固定感染期）。

D. 活动驱动网络上的 SIS 阈值

设活动潜能密度为 $F(a)$，每步活跃节点向 $m$ 个随机节点发出接触，接触持续一步；每次 S–I 接触传播概率为 $\lambda$，感染者每步恢复概率为 $\mu$。记 $I_a$ 为活动性为 $a$ 的类中的感染人数，$N_a=NF(a)$。易感者有两种被感染方式：自己激活并连到感染者（概率 $\approx am\lambda I/N$），或被激活的感染者连到（概率 $\approx m\lambda\int a'I_{a'}da'/N$）。于是 $$I_a^{t+1}-I_a^{t}=-\mu I_a+\lambda m(N_a-I_a)\,a\frac{I}{N}+\lambda m(N_a-I_a)\frac{\int a'I_{a'}\,da'}{N}.$$ 在无病态附近令 $N_a-I_a\approx N_a$，对 $a$ 积分得 $I=\int I_a\,da$ 与 $\Phi=\int aI_a\,da$ 的线性系统： $$\begin{pmatrix}\Delta I\\ \Delta\Phi\end{pmatrix}= \begin{pmatrix}\lambda m\langle a\rangle-\mu & \lambda m\\ \lambda m\langle a^2\rangle & \lambda m\langle a\rangle-\mu\end{pmatrix} \begin{pmatrix}I\\ \Phi\end{pmatrix}.$$ 矩阵的最大特征值为 $\lambda m(\langle a\rangle+\sqrt{\langle a^2\rangle})-\mu$，它大于 0 即为流行条件： $$\frac{\lambda}{\mu}>\frac{1}{m\big(\langle a\rangle+\sqrt{\langle a^2\rangle}\big)}.$$ 若所有人活动性相同（$\langle a^2\rangle=\langle a\rangle^2$），阈值为 $1/(2m\langle a\rangle)$；活动异质性越大，阈值越低。代码示例中，20% 高活跃个体（$a=0.5$）与 80% 低活跃个体（$a=0.05$）的混合使阈值从“同质”估计的 0.354 降到 0.270，模拟的稳态流行率在 0.9 倍理论阈值时为 0、1.2 倍时约为 0.11，转变位置与异质理论一致（代码）。

E. 淬火平均场与非回溯矩阵阈值

QMF：在 $\rho_i=0$ 处线性化个体平均场方程 $\dot\rho_i=-\mu\rho_i+\beta(1-\rho_i)\sum_jA_{ij}\rho_j$，得 $\dot\rho=(\beta A-\mu I)\rho$。$A$ 是对称非负矩阵，最大特征值 $\Lambda_{\max}$ 对应非负特征向量（Perron–Frobenius），故 $\rho$ 沿该方向增长的条件是 $\beta\Lambda_{\max}>\mu$，即 $\lambda_c=1/\Lambda_{\max}$。由 Rayleigh 商，$\Lambda_{\max}\ge\sqrt{k_{\max}}$（取以最大度节点为中心的星形试探向量）且 $\Lambda_{\max}\ge\langle k\rangle$。

非回溯矩阵：消息传递方程在 $u\equiv1$ 附近线性化。令 $u_{i\leftarrow j}=1-\delta_{j\to i}$，保留一阶项得 $\delta_{j\to i}=T\sum_{l\in\partial j\setminus i}\delta_{l\to j}$，即 $\delta=TB\delta$。非零解出现的条件是 $T\rho(B)>1$。

Ihara–Bass 化简：设 $v$ 是 $B$ 的特征向量，特征值为 $\nu\ne0$。由 $B$ 的定义，$\nu v_{j\to i}=\sum_{l\in\partial j}v_{l\to j}-v_{i\to j}$。定义两个节点量：流入之和 $x_j=\sum_{l\in\partial j}v_{l\to j}$，流出之和 $y_j=\sum_{i\in\partial j}v_{j\to i}$。

• 对 $i\in\partial j$ 求和：$\nu y_j=k_jx_j-x_j$，即 $\nu y=(D-I)x$；
• 对固定的 $i$、$j\in\partial i$ 求和：$\nu x_i=\sum_{j\in\partial i}x_j-y_i$，即 $\nu x=Ax-y$。

消去 $y$ 得 $\big(\nu^2I-\nu A+(D-I)\big)x=0$，等价于 $$\begin{pmatrix}A & I-D\\ I & 0\end{pmatrix}\begin{pmatrix}x\\ x/\nu\end{pmatrix}=\nu\begin{pmatrix}x\\ x/\nu\end{pmatrix}.$$ 因此 $\rho(B)$ 可以从这个 $2N\times2N$ 稀疏矩阵的最大特征值求得，计算量与 $A$ 同阶。

F. 目标免疫所需的接种比例

取连续近似 $P(k)=(\gamma-1)k_{\min}^{\gamma-1}k^{-\gamma}$（$k\ge k_{\min}$，$2<\gamma<3$），接种度最高的比例 $g$ 的节点，等价于把度分布在 $k_c$ 处截断： $$g=\int_{k_c}^{\infty}P(k)\,dk=\Big(\frac{k_{\min}}{k_c}\Big)^{\gamma-1}\ \Longrightarrow\ k_c=k_{\min}\,g^{-1/(\gamma-1)}.$$ 忽略“删去枢纽也削减了其邻居的度”这一效应（它会使结果更有利），剩余网络中 $\langle k\rangle\approx k_{\min}\frac{\gamma-1}{\gamma-2}$，$\langle k^2\rangle\approx\frac{\gamma-1}{3-\gamma}k_{\min}^{\gamma-1}k_c^{3-\gamma}$，于是 HMF 阈值 $$\lambda_c(g)=\frac{\langle k\rangle}{\langle k^2\rangle}=\frac{3-\gamma}{(\gamma-2)\,k_{\min}}\;g^{(3-\gamma)/(\gamma-1)} .$$ 令 $\lambda_c(g_c)=\lambda$ 解出 $$g_c=\Big[\frac{(\gamma-2)\,k_{\min}\,\lambda}{3-\gamma}\Big]^{(\gamma-1)/(3-\gamma)} .$$ 算例：$\gamma=2.5$、$k_{\min}=2$，则 $\lambda_c(g)=0.5\,g^{1/3}$；若 $\lambda=0.1$，$g_c=0.2^3=0.008$，即只需接种度最高的约 0.8% 的人。随机免疫则要求 $\lambda(1-g)<\lambda_c(N)$，而 $\lambda_c(N)\to0$，需要几乎全员接种。这一对比与网络鲁棒性、渗流与级联失效中“随机失效鲁棒、蓄意攻击脆弱”是同一现象的两面。

G. 阈值模型的全局级联条件

Watts（2002）考虑配置模型上的阈值模型：节点阈值 $\phi$ 服从分布 $f(\phi)$，节点在已激活邻居比例 $\ge\phi$ 时激活。关键概念是易感节点（vulnerable node）：只要一个邻居激活就会激活的节点，即 $\phi\le1/k$。度为 $k$ 的节点易感的概率为 $\rho_k=\Pr(\phi\le1/k)$。

一个小种子能触发全局级联，当且仅当易感节点构成的子网络存在“渗流”的巨分支。沿一条边到达一个度为 $k$ 的节点的概率为 $kP(k)/\langle k\rangle$，该节点易感的概率为 $\rho_k$，易感时它会把激活继续传给其余 $k-1$ 个邻居。分支过程的平均分支数大于 1，得到级联条件 $$\sum_k k(k-1)\,\rho_k\,P(k)>\langle k\rangle .$$ 算例：所有节点阈值 $\phi=0.18$，则只有 $k\le5$ 的节点易感。对平均度为 $z$ 的泊松网络，条件化为 $z\cdot\Pr(\mathrm{Poisson}(z)\le3)>1$，数值求解得级联窗口约为 $1.02<z<5.76$。$z$ 太小，网络本身不连通；$z$ 太大，几乎没有节点易感（单个邻居的影响 $1/k$ 太小）。窗口上边界附近的级联罕见但一旦发生就覆盖几乎全网，这是“全局级联”难以预测的原因。