03 DFA 算法原理与数学推导

第 3 章 DFA 算法原理与数学推导

本章逐步拆解 DFA 的五步计算流程，并对照 fathon 的 computeFlucVec 参数，建立「公式 ↔ API」的一一映射。

3.1 DFA 算法总览

Detrended Fluctuation Analysis 由 Peng et al. (1994) 提出，Kantelhardt et al. (2002) 系统化。标准流程：

Step 1  构造 profile：Y(k) = Σ_{i=1}^{k} [x(i) - x̄]
Step 2  选定窗口尺度 n，将 profile 分为 N_n 个不重叠（或重叠）段
Step 3  每段内拟合 m 阶多项式 P_ν(i)，计算残差 ε_ν(i) = Y_ν(i) - P_ν(i)
Step 4  计算每段方差并平均，得到 F²(n)，开方得 F(n)
Step 5  对多个 n 重复 Step 2–4，log-log 拟合得 H
​

fathon 中 Step 1 由 fathonUtils.toAggregated 完成；Step 2–5 由 DFA.computeFlucVec + DFA.fitFlucVec 完成。

3.2 Step 1：Profile 构造

给定长度 L 的序列 {x(1), …, x(L)}：

[
\bar{x} = \frac{1}{L}\sum_{t=1}^{L} x(t)
]

[
Y(k) = \sum_{i=1}^{k} [x(i) - \bar{x}], \quad k = 1, \ldots, L
]

代码对应：

from fathon import fathonUtils as fu
import numpy as np

x = np.random.randn(5000)
Y = fu.toAggregated(x)  # 等价于 np.cumsum(x - np.mean(x))
​

subtractMean 仅去均值不累积；toAggregated = 去均值 + 累积和。

3.3 Step 2：窗口分割

对窗口大小 n（n ≪ L），将 profile 分为：

[
N_n = \lfloor L / n \rfloor
]

个不重叠段（默认）。每段 ν 包含样本 Y_ν(1), …, Y_ν(n)。

反向分段（revSeg=True）：从序列末尾再分一遍，共 2N_n 段，波动取平均。这减少端点效应——序列开头/结尾不足一整窗的部分被更好利用。

fathon 参数映射：

3.4 Step 3：局部去趋势

在第 ν 段内，拟合 polOrd 阶多项式：

[
P_\nu(i) = a_0 + a_1 i + a_2 i^2 + \cdots + a_m i^m, \quad m = \texttt{polOrd}
]

最小二乘拟合后，残差：

[
\epsilon_\nu(i) = Y_\nu(i) - P_\nu(i)
]

polOrd 选择指南（详见第 9 章）：

3.5 Step 4：波动函数 F(n)

每段的去趋势方差：

[
F^2_\nu(n) = \frac{1}{n}\sum_{i=1}^{n} \epsilon_\nu^2(i)
]

所有段平均后开方：

[
F(n) = \sqrt{\frac{1}{2N_n}\sum_{\nu=1}^{2N_n} F^2_\nu(n)} \quad \text{（revSeg 时）}
]

或仅正向 N_n 段（revSeg=False）。

标度假设：若序列具有长记忆标度，则：

[
F(n) \sim n^H
]

3.6 Step 5：双对数拟合求 H

对窗口数组 {n_j} 及对应 {F(n_j)}，在 log-log 空间线性回归：

[
\log F(n_j) = H \cdot \log n_j + c
]

斜率 H 即 Hurst 指数。fathon 的 fitFlucVec 实现此拟合，默认底数为 e：

H, intercept = dfa.fitFlucVec(logBase=np.e)
# log(F) = H * ln(n) + intercept  （自然对数）
​

若改用 logBase=2，斜率 H 不变（对数底变换只影响截距）。

3.7 无偏 DFA（unbiased=True）

v1.3.2 引入无偏版本（Alvarez-Ramirez et al., 2018），适用于短序列：

n, F = dfa.computeFlucVec(wins, unbiased=True)
# 注意：unbiased=True 时 revSeg 被忽略
​

无偏 DFA 修正有限样本偏差，使短序列 H 估计更接近渐近值。序列长度 L < 1000 时建议尝试。

3.8 手算迷你示例（理解 F(n) 含义）

下面用纯 NumPy 实现单窗口 n=100 的核心逻辑，帮助对照 fathon 黑盒：

import numpy as np

def manual_F(profile, n, polOrd=1):
    """单尺度 n 的 DFA 波动（简化版，仅正向分段）"""
    L = len(profile)
    N_n = L // n
    F2_list = []
    for nu in range(N_n):
        seg = profile[nu * n : (nu + 1) * n]
        t = np.arange(n, dtype=float)
        # 多项式拟合：最小二乘
        coeffs = np.polyfit(t, seg, polOrd)
        trend = np.polyval(coeffs, t)
        residual = seg - trend
        F2_list.append(np.mean(residual ** 2))
    return np.sqrt(np.mean(F2_list))

np.random.seed(0)
x = np.random.randn(4096)
Y = np.cumsum(x - np.mean(x))

for n in [16, 64, 256, 1024]:
    F = manual_F(Y, n)
    print(f"n={n:5d}  F(n)={F:.4f}")
​

运行后 log F 对 log n 近似线性，斜率 ≈ 0.5（白噪声）。

3.9 完整 fathon 对照实验

import numpy as np
import fathon
from fathon import fathonUtils as fu

np.random.seed(0)
L = 4096
profile = fu.toAggregated(np.random.randn(L))
wins = np.array([16, 64, 256, 1024])

dfa = fathon.DFA(profile)
n, F = dfa.computeFlucVec(wins, revSeg=False, polOrd=1)
H, _ = dfa.fitFlucVec()

print("窗口 n:", n)
print("F(n):  ", np.round(F, 4))
print("H:     ", round(H, 4))

# 手算验证斜率
log_n = np.log(n)
log_F = np.log(F)
H_manual = np.polyfit(log_n, log_F, 1)[0]
print("手算斜率:", round(H_manual, 4))
​

fathon 的 H 与 np.polyfit 手算斜率应高度一致。

3.10 计算复杂度与 fathon 实现

fathon 核心用 Cython + C + GSL 实现，并支持 OpenMP 多线程（含 Windows）。

长序列 + 多 q 的 MFDFA 建议先用较大步长探索 wins，再细化（第 8 章）。

3.11 常见误区

3.12 本章小结

x(t) ──toAggregated──► Y(t) profile
                          │
              computeFlucVec(winSizes, polOrd, revSeg)
                          │
                    F(n) for each n
                          │
                   fitFlucVec() ──► H, intercept
​

下一章聚焦 DFA 类的完整 API 与可视化模板。

3.13 动手实验

1. 运行 3.8 节手算代码，与 3.9 节 fathon 输出对比 F(n) 数值。
2. 设 revSeg=True/False 各跑一次，记录 H 差异。
3. 对长度 L=500 的序列分别用 unbiased=False/True，比较 H（短序列实验）。

3.14 配套可运行示例

本章的「手算 vs 加速版对照」已封装成示例，对照阅读效果最佳：

💡 强烈建议运行一遍「纯 NumPy 手算 DFA」示例：自己实现一遍 DFA 五步，再与加速版结果对比，你会真正理解 F(n) 不是黑盒。

下一章聚焦 DFA 类的完整 API 与可视化模板。