数学推导与算法

本节推导滑动窗口相关的方差与窗口长度的关系，给出 Allen 等（2014）连接状态流程、相位随机化替代数据、HMM、多层模块度和 ICC 的要点。符号：$N$ 个脑区，$T$ 个时间点，窗长 $w$（时间点数），步长 $s$，第 $m$ 个窗口的相关矩阵为 $R^{(m)}$。

A. 滑动窗口估计量：方差从哪里来

在窗口 $m$ 内，脑区 $i,j$ 的估计量是样本相关 $r_{ij}^{(m)}$。先假设真实相关恒为 $\rho$（完全平稳、没有任何动态）。

独立样本情形：由 Fisher 变换的方差稳定化结果（推导见脑连接组的数学推导），$z^{(m)}=\operatorname{artanh}(r^{(m)})$ 近似服从 $\mathcal N\big(\operatorname{artanh}\rho,\ 1/(w-3)\big)$。因此跨窗口的波动方差

$$\operatorname{Var}_m\big(z^{(m)}\big)\approx\frac{1}{w-3}$$

完全由窗长决定，与大脑是否“动态”无关。窗长减半，方差约加倍。

自相关情形：BOLD 信号相邻时间点高度相关。设两信号的自相关函数分别为 $\rho_{ii}(\ell)$、$\rho_{jj}(\ell)$（$\ell$ 为滞后）。在 $\rho=0$ 的零假设下，Bartlett 公式给出

$$\operatorname{Var}(r)\approx\frac{1}{w}\sum_{\ell=-\infty}^{\infty}\rho_{ii}(\ell)\,\rho_{jj}(\ell)\;\equiv\;\frac{1}{w_{\text{eff}}}.$$

若两信号都是系数为 $\phi$ 的 AR(1) 过程，$\rho(\ell)=\phi^{|\ell|}$，则

$$\sum_{\ell}\phi^{2|\ell|}=1+\frac{2\phi^2}{1-\phi^2}=\frac{1+\phi^2}{1-\phi^2},\qquad w_{\text{eff}}=w\cdot\frac{1-\phi^2}{1+\phi^2}.$$

例如 $\phi=0.7$ 时 $w_{\text{eff}}\approx0.34\,w$：一个 30 点的窗口只相当于约 10 个独立样本，$z$ 的标准差从 0.19 增大到约 0.3 以上。

偏差一侧：短窗有方差问题，长窗则有偏差问题。把两条信号在窗内标准化后，窗口相关近似等于两信号逐点乘积 $x_{it}x_{jt}$（即边时间序列）的滑动平均。长度为 $w$ 的滑动平均是一个低通滤波器，第一个零点在频率 $1/w$ 处，所以快于约 $1/w$ 的真实连接变化会被抹平。窗长同时控制“噪声有多大”和“能看到多快的变化”，这是典型的偏差–方差权衡；Leonardi & Van De Ville（2015）的 $w\ge1/f_{\min}$ 建议，则是为避免信号本身的低频成分漏入相关估计而给出的下限。

B. 连接状态流程（Allen 等, 2014）

算法 1  基于 k-means 的动态连接状态
输入：S 名被试的节点时间序列 x_s ∈ R^{N×T}；窗长 w；步长 s；候选状态数 K_set
输出：K 个状态中心；每名被试的状态序列及占时比例、停留时间、转移矩阵
1. 对每名被试 s、每个窗口 m：
2.     用锥形窗（矩形窗与高斯核卷积）加权，计算 N×N 连接矩阵 R_s^(m)
       （原文用 L1 正则化的逆协方差估计，以缓解窗内样本少的问题）
3.     Fisher z 变换，取上三角拉直为向量 v_s^(m)，长度 N(N−1)/2
4. 选“代表窗口”：每名被试取连接方差序列的局部极大处的窗口
5. 对代表窗口做 k-means（原文用 L1 距离），用肘部准则在 K_set 中选 K
6. 以第 5 步中心为初值，对全部窗口做 k-means，得到最终中心 c_1..c_K
7. 把每名被试的窗口按最终标签排成状态序列 ℓ_s(m)
8. 由 ℓ_s(m) 计算：占时比例、平均停留时间、K×K 转移概率矩阵
复杂度：建窗 O(S·M·N²·w)，M ≈ (T−w)/s+1；k-means 每次迭代 O(S·M·K·N²)

关键参数：$w$（原文 22 个 TR，约 44 秒）、$K$（原文选出 7）、距离度量与初始化。三者都应做敏感性分析。第 5 步的肘部准则只是经验规则；可用不同随机初值重复聚类，检查状态中心的匹配程度（见下面 F）。

C. 相位随机化替代数据

要检验“波动是否超出平稳预期”，需要一个保留功率谱和静态协方差、但破坏非平稳结构的零模型。多元情形的关键是所有通道加同一组随机相位：离散傅里叶变换后，每个频率 $f$ 处的交叉谱 $X_i(f)X_j^*(f)$ 在 $X_i,X_j$ 同乘 $e^{i\varphi_f}$ 后不变，所以每个信号的功率谱和信号间的线性相关结构（交叉谱）都保留下来。

算法 2  多元相位随机化替代数据（Prichard & Theiler 思路）
输入：X ∈ R^{N×T}（各通道已去均值）；替代数 B（如 500–1000）；检验统计量 g(·)
输出：g 的零分布及经验 p 值
1. g_obs ← g(X)，例如“窗口相关时间序列的方差”或状态切换次数
2. F ← FFT(X)，按行（每通道）计算
3. 对 b = 1..B：
4.     为正频率 f = 1..⌊(T−1)/2⌋ 抽相位 φ_f ~ U(0, 2π)
5.     负频率取 −φ_f，保证逆变换为实数；直流（及偶数 T 的奈奎斯特频率）相位保持不变
6.     F_b(:, f) ← F(:, f) · exp(i φ_f)     （所有通道共用同一个 φ_f）
7.     X_b ← real(IFFT(F_b))；g_b ← g(X_b)
8. p ← (1 + #{g_b ≥ g_obs}) / (B + 1)
复杂度：O(B · N · T log T)，再加上每次计算 g 的开销（通常更大）

D. 边时间序列恒等式

设 $\bar x_i=\frac1T\sum_tx_i(t)$，$\sigma_i^2=\frac1T\sum_t(x_i(t)-\bar x_i)^2$。皮尔逊相关

$$r_{ij}=\frac{\frac1T\sum_t(x_i-\bar x_i)(x_j-\bar x_j)}{\sigma_i\sigma_j}=\frac1T\sum_t z_i(t)z_j(t)=\frac1T\sum_te_{ij}(t).$$

（若用 $T-1$ 作分母的样本标准差，只差一个常数因子 $(T-1)/T$。）由此可以理解静态零模型的论证：若 $\mathbf z_t$ 是独立同分布的高斯向量，协方差为 $R$，则 $e_{ij}(t)$ 的均值为 $r_{ij}$，方差为 $1+r_{ij}^2$（高斯四阶矩的 Isserlis 定理）；全部边的平方和 $\sum_{i<j}e_{ij}(t)^2$ 较大的帧，也就是 $\mathbf z_t$ 沿 $R$ 的主要特征方向取值较大的帧，它们的外积 $\mathbf z_t\mathbf z_t^\top$ 自然更接近 $R$ 的主要结构。所以“高幅度帧更像静态 FC”在完全平稳的数据中也会出现。

E. HMM 的前向–后向与 EM

记发射概率 $b_t(k)=p(\mathbf x_t\mid s_t=k)$。

• 前向：$\alpha_1(k)=\pi_kb_1(k)$，$\alpha_t(k)=b_t(k)\sum_j\alpha_{t-1}(j)A_{jk}$。实际计算中每步归一化：$c_t=\sum_k\alpha_t(k)$，$\hat\alpha_t=\alpha_t/c_t$，对数似然为 $\sum_t\log c_t$。
• 后向：$\hat\beta_T(k)=1$，$\hat\beta_t(k)=\frac{1}{c_{t+1}}\sum_lA_{kl}b_{t+1}(l)\hat\beta_{t+1}(l)$。
• 后验：$\gamma_t(k)=\hat\alpha_t(k)\hat\beta_t(k)$；$\xi_t(k,l)=\hat\alpha_t(k)A_{kl}b_{t+1}(l)\hat\beta_{t+1}(l)/c_{t+1}$。
• M 步：$\pi_k\leftarrow\gamma_1(k)$；$A_{kl}\leftarrow\sum_t\xi_t(k,l)/\sum_t\sum_{l'}\xi_t(k,l')$；$\Sigma_k\leftarrow\sum_t\gamma_t(k)\mathbf x_t\mathbf x_t^\top/\sum_t\gamma_t(k)$（零均值模型）。

每次迭代复杂度 $O(TK^2+TKN^2)$（另加求逆和行列式的 $O(KN^3)$）。EM 只保证收敛到局部极大，因此需要多次初始化并比较似然；状态数的选择见状态数的选择与稳定性。

F. 多层模块度、局部移动与灵活性

把多层模块度写成超节点 $(i,l)$ 上的二次型：定义超模块度矩阵

$$\mathcal B_{il,jr}=\Big(A_{ijl}-\gamma_l\frac{k_{il}k_{jl}}{2m_l}\Big)\delta_{lr}+\delta_{ij}\,\omega_{jlr},$$

则 $Q_{\text{ml}}=\frac{1}{2\mu}\sum\mathcal B_{il,jr}\,\delta(g_{il},g_{jr})$，形式与单层模块度 $Q=\frac{1}{2m}\sum B_{ij}\delta(g_i,g_j)$ 完全相同。因此任何以模块度矩阵为输入的算法都可直接使用。

算法 4  超模块度矩阵上的 Louvain（GenLouvain 思路）
输入：NL×NL 矩阵 𝓑；输出：每个超节点的社团标签
1. 每个超节点自成一个社团
2. 局部移动：随机顺序访问节点 v，计算把 v 移入各社团 c 后的得分 Σ_{u∈c, u≠v} 𝓑_{vu}，
   移到得分最高的社团；重复直到没有节点移动
3. 聚合：把每个社团合并为一个新节点，新矩阵 𝓑' = Pᵀ 𝓑 P（P 为成员指示矩阵）
4. 在 𝓑' 上重复 2–3，直到不再合并
5. 多次随机运行，比较 Q 并求共识划分

灵活性 $f_i$ 与模块忠诚度 $P_{ij}$ 都可以从最终标签 $g_{il}$ 直接计算（定义见网络灵活性与动态社团指标）。

G. 状态匹配（匈牙利算法）

设两次聚类的中心为 $\{\mathbf c_k\}$ 与 $\{\mathbf c'_l\}$，相似度矩阵 $S_{kl}=\operatorname{corr}(\mathbf c_k,\mathbf c'_l)$。求置换 $\sigma$ 使 $\sum_kS_{k\sigma(k)}$ 最大，这是一个线性指派问题，匈牙利算法在 $O(K^3)$ 时间内求出最优解。报告 $\min_kS_{k\sigma(k)}$（最差匹配）比报告平均值更能暴露不稳定的状态。

H. ICC 的方差分析形式

把 $n\times k$ 数据矩阵 $Y$（被试×次数）的总平方和分解为被试间 $SS_R$、次数间 $SS_C$ 与残差 $SS_E$：$MS_R=SS_R/(n-1)$，$MS_C=SS_C/(k-1)$，$MS_E=SS_E/[(n-1)(k-1)]$。则

$$\mathrm{ICC}(3,1)=\frac{MS_R-MS_E}{MS_R+(k-1)MS_E},\qquad \mathrm{ICC}(2,1)=\frac{MS_R-MS_E}{MS_R+(k-1)MS_E+\frac{k}{n}(MS_C-MS_E)}.$$

当测量误差方差 $\sigma_e^2$ 远大于被试间真实方差 $\sigma_b^2$ 时，$\mathrm{ICC}\approx\sigma_b^2/(\sigma_b^2+\sigma_e^2)\to0$：这正是短窗口动态指标信度低的原因。