10 进阶拓展方向

文档摘要

第 10 章进阶拓展方向掌握 fathon 五类算法后，本章指向合成数据验证、领域应用、与其他实现对比及研究前沿。 10.1 fBm 合成数据验证验证 DFA 实现的金标准是分数布朗运动（fBm），其理论 Hurst 已知。Davies-Harte 算法可在 O(n log n) 生成 fBm 样本：注意：完整 Davies-Harte 实现较复杂；生产验证可用 R 的或 MATLAB 工具箱对照。上例 Cholesky 版仅适用于中等 n。 10.2 金融时间序列工作流 10.

第 10 章进阶拓展方向

掌握 fathon 五类算法后，本章指向合成数据验证、领域应用、与其他实现对比及研究前沿。

10.1 fBm 合成数据验证

验证 DFA 实现的金标准是 分数布朗运动（fBm），其理论 Hurst 已知。Davies-Harte 算法可在 O(n log n) 生成 fBm 样本：


import numpy as np

def davies_harte_fbm(n, H, seed=None):
    """
    生成 fBm 在时刻 1..n 的值（简化 Davies-Harte）。
    n: 样本点数（建议 2 的幂次）
    H: Hurst 指数
    """
    rng = np.random.default_rng(seed)
    if n % 2 != 0:
        n = n + 1
    m = n // 2
    # 频域权重
    k = np.arange(1, m + 1)
    fk = np.zeros(m + 1, dtype=complex)
    sigma = np.sqrt(0.5 * np.array([
        0.0,
        *[0.5 * (k_i ** (-2 * H - 1)) for k_i in k],
    ]))
    # 构造周期延拓的 circulant 嵌入（简化版）
    g = rng.standard_normal(m + 1) + 1j * rng.standard_normal(m + 1)
    fk = sigma * g
    fk[0] = 0
    # 逆 FFT 得 fGn 相关结构，再积分得 fBm 增量
    # 此处用 Cholesky 近似（小 n 演示）
    t = np.arange(1, n + 1)
    cov = 0.5 * (np.add.outer(t, t) ** (2 * H) -
                 np.abs(np.subtract.outer(t, t)) ** (2 * H))
    cov += 1e-10 * np.eye(n)
    L_chol = np.linalg.cholesky(cov)
    fbm = L_chol @ rng.standard_normal(n)
    return fbm

# 验证 pipeline
import fathon
from fathon import fathonUtils as fu

for H_true in [0.3, 0.5, 0.7, 0.9]:
    fbm = davies_harte_fbm(4096, H_true, seed=42)
    increments = np.diff(np.concatenate([[0], fbm]))
    profile = fu.toAggregated(increments)
    wins = fu.linRangeByStep(10, 1024, step=8)
    dfa = fathon.DFA(profile)
    n, F = dfa.computeFlucVec(wins, revSeg=True, polOrd=1)
    H_est, _ = dfa.fitFlucVec()
    print(f"H_true={H_true:.1f}  H_est={H_est:.3f}  err={abs(H_est-H_true):.3f}")

注意：完整 Davies-Harte 实现较复杂；生产验证可用 R 的 somefbm 或 MATLAB 工具箱对照。上例 Cholesky 版仅适用于中等 n。

10.2 金融时间序列工作流


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

# 模拟日对数收益率
np.random.seed(0)
T = 252 * 10  # 10 年日频
returns = np.random.randn(T) * 0.01
returns[:T//2] *= 1.5  # 前半波动更大（多重分形 hint）

profile = fu.toAggregated(returns)
wins = fu.powRangeByCount(4, int(np.log2(T // 4)), count=40, base=2).astype(int)

# 1. 全局 DFA
dfa = fathon.DFA(profile)
n, F = dfa.computeFlucVec(wins, revSeg=True)
H, _ = dfa.fitFlucVec()
print(f"收益率 DFA H = {H:.4f}")

# 2. MFDFA 检测波动率 regime
q_list = np.arange(-4, 5, 0.5)
mfdfa = fathon.MFDFA(profile)
n2, F2 = mfdfa.computeFlucVec(wins, q_list, revSeg=True)
h_q, _ = mfdfa.fitFlucVec()
alpha, f_alpha = mfdfa.computeMultifractalSpectrum()
print(f"Δh = {h_q.max()-h_q.min():.3f}, Δα = {alpha.max()-alpha.min():.3f}")

# 3. 两资产交叉 DCCA（演示用第二列随机）
r2 = returns * 0.6 + np.random.randn(T) * 0.008
p1, p2 = fu.toAggregated(returns), fu.toAggregated(r2)
dcca = fathon.DCCA(p1, p2)
n3, F3 = dcca.computeFlucVec(wins[:20], revSeg=True)
H_x, _ = dcca.fitFlucVec()
print(f"交叉 DCCA H = {H_x:.4f}")

10.3 生理信号：RR 间期

心率 RR 间期序列常呈长记忆与多分形：

步骤	操作
预处理	去除异位搏、插值缺失
增量	RR 一阶差分或直接 toAggregated(RR-mean(RR))
DFA	polOrd=1～2，关注 H 是否 > 0.5
MFDFA	健康 vs 病理的 Δα 对比
HT	检测自主神经调节的时变 persistence

采样率与长度要求：通常 L > 2000 搏。

10.4 fathon vs 其他 Python/R 实现

库	语言	特点
fathon	Python	DFA+MFDFA+DCCA+MFDCCA+HT 全家桶，Cython 加速
MFDFA (LRydin)	Python	专注 MFDFA，EMD 选项，向量化多项式
nolds	Python	多种混沌/非线性指标含 DFA
R fractal	R	经典 fractal 分析

选型建议：

需要 DCCA/MFDCCA/HT → fathon 首选
只需 MFDFA + EMD → 可考虑 MFDFA 包
快速 H 估计 嵌入 ML pipeline → fathon DFA 或 nolds

10.5 与机器学习集成

DFA/MFDFA 特征可作为时序分类器的手工特征：


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

def extract_features(increments, q_list=None):
    """从增量序列提取标度特征向量"""
    L = len(increments)
    profile = fu.toAggregated(increments)
    wins = fu.linRangeByStep(10, L // 4, step=max(1, L // 200))

    feats = {}

    dfa = fathon.DFA(profile)
    n, F = dfa.computeFlucVec(wins, revSeg=True, polOrd=1)
    H, _ = dfa.fitFlucVec()
    feats['H'] = H

    if q_list is not None and L >= 2000:
        mfdfa = fathon.MFDFA(profile)
        _, Fq = mfdfa.computeFlucVec(wins, q_list, revSeg=True)
        hq, _ = mfdfa.fitFlucVec()
        feats['h_q_mean'] = np.mean(hq)
        feats['h_q_std'] = np.std(hq)
        feats['delta_h'] = hq.max() - hq.min()

    return feats

# 演示
inc = np.random.randn(3000)
q = np.arange(-2, 3, 1.0)
print(extract_features(inc, q_list=q))

可将 H、delta_h 等与 CNN/LSTM 特征拼接，或用于 XGBoost 等 tabular 模型。

10.6 研究前沿与扩展阅读

方向	关键词
高阶 DFA	三阶矩、不对称波动
多元 DFA	多通道联合标度
小波领导	Wavelet Leaders 多分形
自适应去趋势	EMD/SSA 替代多项式
非平稳 MFDFA	滑动 MFDFA

原始论文入口：

Peng et al. (1994) — DFA 起源
Kantelhardt et al. (2002) — MFDFA
Podobnik & Stanley (2008) — DCCA
Ihlen (2012) — 时变 Hurst
Bianchi (2020) — fathon JOSS 论文

10.7 贡献与反馈

fathon 开源（GPL-3.0）。发现问题可至 GitHub 仓库提交 Issue。引用请用：

Bianchi, S. (2020). fathon: A Python package for a fast computation of detrended fluctuation analysis and related algorithms. Journal of Open Source Software, 5(45), 1828.

10.8 全教程知识地图


第 0～2 章   概念与动机
第 3～4 章   DFA 原理与 API  ← 核心基础
第 5 章      MFDFA 多分形
第 6 章      DCCA / MFDCCA
第 7 章      HT 时变 Hurst
第 8 章      工具与窗口
第 9 章      工程实践
第 10 章     进阶（本章）
附录 A～C    术语、速查、排错

10.9 动手实验（毕业项目建议）

合成验证：用 10.1 节对不同 H 的 fBm 估计误差表。
公开数据集：下载一组金融/生理公开数据，完成第 9 章报告模板。
方法对比：同一序列分别用 fathon 与 nolds（若安装）估计 H，对比差异。
分类小项目：两类序列（如 persistent vs white），用 MFDFA 的 Δh 做简单分类器。

10.10 配套可运行示例

「分式布朗运动合成与标定」示例用 Davies-Harte 法生成已知 Hurst 的 fBm，再用 DFA 反推 H，验证估计准确度。下面是完整脚本：


"""分式布朗运动（fBm）合成与 DFA 标定完整示例。

用 Davies-Harte 近似方法生成给定 Hurst 的 fBm，然后用 DFA 反推 H，
验证估计流程在已知答案的合成数据上的准确度。
依赖：pip install fathon numpy
"""
import numpy as np
import fathon
from fathon import fathonUtils as fu


def fbm_davies_harte(H, n, seed=None):
    """分式布朗运动的 Davies-Harte 近似生成。

    fBm 自协方差：γ(k) = 0.5*(|k+1|^(2H) - 2|k|^(2H) + |k-1|^(2H))
    利用循环嵌入 + FFT 得到长度为 n 的 fBm 样本。
    """
    rng = np.random.default_rng(seed)
    k = np.arange(n)
    gamma = 0.5 * (np.abs(k + 1) ** (2 * H)
                   - 2 * np.abs(k) ** (2 * H)
                   + np.abs(k - 1) ** (2 * H))
    circ = np.concatenate([gamma, gamma[-2:0:-1]])
    spectrum = np.fft.fft(circ).real
    spectrum = np.maximum(spectrum, 0)
    sqrt_spectrum = np.sqrt(spectrum)
    complex_normal = (rng.standard_normal(len(circ))
                      + 1j * rng.standard_normal(len(circ))) / np.sqrt(2)
    freq = sqrt_spectrum * complex_normal
    fbm_sample = np.fft.ifft(freq).real[:n]
    return fbm_sample / np.std(fbm_sample) * np.sqrt(gamma[0])


def main():
    wins = fu.linRangeByStep(10, 1500, step=10)
    print("=== 用 DFA 反推已知 H 的 fBm ===")
    print(f"{'设定H':>8} {'估计H':>8} {'误差':>8}")
    for H_true in [0.3, 0.5, 0.7, 0.9]:
        bm = fbm_davies_harte(H_true, 4000, seed=int(H_true * 100))
        profile = fu.toAggregated(bm)
        dfa = fathon.DFA(profile)
        dfa.computeFlucVec(wins, revSeg=True)
        H_est, _ = dfa.fitFlucVec()
        err = H_est - H_true
        print(f"{H_true:8.2f} {H_est:8.4f} {err:+8.4f}")


if __name__ == "__main__":
    main()

💡 这是检验你整个 DFA 流程（profile / 窗口 / polOrd）是否正确最直接的方法——如果连已知答案都估不对，真实数据上更不能信。若偏差大于 0.05，回头检查窗口上限、polOrd、是否对增量而非原序列做 toAggregated。

10.11 结语

DFA 及其扩展是理解复杂时间序列长记忆与多分形结构的利器。fathon 将五种算法统一在同一 API 风格下，配合本教程的「原理 → 代码 → 工程」路径，你应能：

对任意一维增量序列估计 Hurst 与 h(q) 谱
分析双序列交叉标度与显著性
在合适的场景使用时变 HT
避开窗口、polOrd、profile 等常见坑

祝你在灏天文库发布顺利，探索愉快。