05 MFDFA 多分形扩展

第 5 章 MFDFA 多分形扩展

当单一 Hurst 指数 H 无法刻画序列在不同波动幅度上的标度差异时，需要 MFDFA（Multifractal Detrended Fluctuation Analysis）。本章讲解 h(q)、τ(q)、f(α) 的含义与 fathon 完整用法。

5.1 从 DFA 到 MFDFA

DFA 的波动函数 F(n) 描述整体标度。MFDFA（Kantelhardt et al., 2002）引入 q 阶广义波动函数：

[
F_q(n) = \left{ \frac{1}{2N_n} \sum_{\nu=1}^{2N_n} [F^2_\nu(n)]^{q/2} \right}^{1/q}
]

广义 Hurst 指数 h(q) 由标度律定义：

[
F_q(n) \sim n^{h(q)}
]

• 若 h(q) 为常数（≈ H）：单标度（monofractal）
• 若 h(q) 随 q 变化：多分形（multifractal）

5.2 质量指数 τ(q) 与多重分形谱 f(α)

质量指数：

[
\tau(q) = q \cdot h(q) - 1
]

奇异性强度与谱（Legendre 变换）：

[
\alpha = h(q) + q \frac{dh(q)}{dq}
]

[
f(\alpha) = q[\alpha - h(q)] + 1
]

f(α) 谱的宽度 Δα = α_max - α_min 衡量多分形强度：越宽，多分形性越强。

fathon 提供：

• computeMassExponents() → τ(q)
• computeMultifractalSpectrum() → α, f(α)

5.3 fathon.MFDFA API 概览

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

profile = fu.toAggregated(np.random.randn(10000))
mfdfa = fathon.MFDFA(profile)

wins = fu.linRangeByStep(10, 2000)
q_list = np.arange(-3, 4, 0.5)   # q 从 -3 到 3.5，步长 0.5

n, F = mfdfa.computeFlucVec(wins, q_list, polOrd=1, revSeg=True)
h_q, h_intercept = mfdfa.fitFlucVec()
tau = mfdfa.computeMassExponents()
alpha, mfSpect = mfdfa.computeMultifractalSpectrum()
​

方法一览

5.4 完整可运行示例 + 四联图

import numpy as np
import fathon
from fathon import fathonUtils as fu
import matplotlib.pyplot as plt

np.random.seed(42)
L = 15000

# 合成「多分形」风格序列：对数正态乘性级联的简化近似
# 用不同方差的分段白噪声拼接，使 h(q) 随 q 变化
segments = []
for scale in [0.5, 1.0, 2.0, 4.0]:
    segments.append(np.random.randn(L // 4) * scale)
increments = np.concatenate(segments)
profile = fu.toAggregated(increments)

wins = fu.linRangeByStep(10, L // 4, step=5)
q_list = np.arange(-4, 5, 0.5)

mfdfa = fathon.MFDFA(profile)
n, F = mfdfa.computeFlucVec(wins, q_list, revSeg=True, polOrd=1)
h_q, _ = mfdfa.fitFlucVec()
tau = mfdfa.computeMassExponents()
alpha, mfSpect = mfdfa.computeMultifractalSpectrum()

# ── 四联图 ──
fig, axes = plt.subplots(2, 2, figsize=(11, 9))

# 1. F_q(n) 标度（选几个 q）
for qi in [0, 2, 4, -2]:
    idx = np.argmin(np.abs(q_list - qi))
    axes[0, 0].loglog(n, F[idx], 'o-', ms=2, label=f'q={q_list[idx]:.1f}')
axes[0, 0].set_xlabel('n'); axes[0, 0].set_ylabel('F_q(n)')
axes[0, 0].set_title('广义波动函数'); axes[0, 0].legend(fontsize=8)

# 2. h(q)
axes[0, 1].plot(q_list, h_q, 'b-o', ms=4)
axes[0, 1].axhline(0.5, color='gray', ls='--', alpha=0.5)
axes[0, 1].set_xlabel('q'); axes[0, 1].set_ylabel('h(q)')
axes[0, 1].set_title('广义 Hurst 指数')

# 3. τ(q)
axes[1, 0].plot(q_list, tau, 'g-s', ms=4)
axes[1, 0].set_xlabel('q'); axes[1, 0].set_ylabel('τ(q)')
axes[1, 0].set_title('质量指数')

# 4. f(α) 谱
axes[1, 1].plot(alpha, mfSpect, 'r-^', ms=4)
axes[1, 1].set_xlabel('α'); axes[1, 1].set_ylabel('f(α)')
axes[1, 1].set_title(f'多重分形谱 (Δα={alpha.max()-alpha.min():.3f})')

plt.tight_layout()
plt.show()

print(f"h(2) ≈ {h_q[np.argmin(np.abs(q_list-2))]:.3f}  (应对标 DFA 的 H)")
print(f"Δh = {h_q.max()-h_q.min():.3f}")
​

5.5 q 列表设计原则

计算成本：F 的 shape 为 (len(q), len(wins))，q 和 wins 都加倍则内存与 CPU 近似乘积增长。

5.6 单标度 vs 多分形 判据

5.7 MFDFA 与 DFA 结果对照

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

np.random.seed(0)
profile = fu.toAggregated(np.random.randn(8000))
wins = fu.linRangeByStep(10, 2000)

# DFA
dfa = fathon.DFA(profile)
_, F_dfa = dfa.computeFlucVec(wins, revSeg=True)
H_dfa, _ = dfa.fitFlucVec()

# MFDFA q=2 应对齐
mfdfa = fathon.MFDFA(profile)
q2_idx = 0  # 若 q_list 第一个为 2
q_list = np.array([2.0])
_, F_mf = mfdfa.computeFlucVec(wins, q_list, revSeg=True)
h2, _ = mfdfa.fitFlucVec()

print(f"DFA  H = {H_dfa:.4f}")
print(f"MFDFA h(2) = {h2[0]:.4f}")
print(f"差异 = {abs(H_dfa - h2[0]):.4f}")
​

两者应非常接近（数值误差 < 0.01 量级）。

5.8 典型报错场景

5.9 本章小结

profile ──► MFDFA ──► computeFlucVec(wins, qList)
                          │
                    fitFlucVec ──► h(q)
                          │
              computeMassExponents ──► τ(q)
                          │
         computeMultifractalSpectrum ──► α, f(α)
​

MFDFA 是 DFA 的自然推广：q=2 退化为 DFA；h(q) 谱宽揭示多分形性。

5.10 动手实验

1. 对白噪声跑 MFDFA，验证 h(q) 是否近似常数 0.5。
2. 将 q 步长从 0.5 改为 0.1，观察 h(q) 平滑度与运行时间。
3. 比较第 5.4 节「拼接序列」与纯白噪声的 Δα。

5.11 配套可运行示例

本章四联图与 DFA/MFDFA 对照封装在「MFDFA 多分形四联图」示例中：同时跑白噪声（单标度对照）与拼接序列（近似多分形），输出 h(q)、τ(q)、f(α) 三套指标与 Δh、Δα。下面是可直接运行的完整脚本：

"""MFDFA 多分形分析完整示例。

对白噪声（期望单标度）与拼接序列（近似多分形）跑 MFDFA，
输出广义 Hurst h(q)、质量指数 τ(q)、多重分形谱 f(α)。
依赖：pip install fathon numpy
"""
import numpy as np
import fathon
from fathon import fathonUtils as fu


def make_multifractal_like(n, seed=0):
    """合成近似多分形序列：不同尺度分段白噪声拼接，使 h(q) 随 q 变化。"""
    rng = np.random.default_rng(seed)
    parts = []
    for scale in [0.5, 1.0, 2.0, 4.0]:
        parts.append(rng.standard_normal(n // 4) * scale)
    return np.concatenate(parts)


def run_mfdfa(profile, wins, q_list, title):
    mfdfa = fathon.MFDFA(profile)
    n, F = mfdfa.computeFlucVec(wins, q_list, revSeg=True, polOrd=1)
    h_q, _ = mfdfa.fitFlucVec()
    tau = mfdfa.computeMassExponents()
    alpha, mfSpect = mfdfa.computeMultifractalSpectrum()
    h2 = h_q[np.argmin(np.abs(q_list - 2))]
    dh = h_q.max() - h_q.min()
    print(f"\n[{title}]")
    print(f"  h(2) = {h2:.4f}  (应对标 DFA 的 H)")
    print(f"  Δh   = {dh:.4f}")
    print(f"  Δα   = {alpha.max() - alpha.min():.4f}")
    # 判据：Δh < 0.1 视为单标度；> 0.1 视为多分形
    verdict = "多分形" if dh > 0.1 else "单标度"
    print(f"  判定 : {verdict}")


def main():
    wins = fu.linRangeByStep(10, 2000, step=5)
    q_list = np.arange(-4, 5, 0.5)

    # 序列 1：白噪声（单标度对照）
    np.random.seed(7)
    profile_white = fu.toAggregated(np.random.randn(10000))

    # 序列 2：拼接序列（近似多分形）
    inc = make_multifractal_like(15000, seed=11)
    profile_mf = fu.toAggregated(inc)

    run_mfdfa(profile_white, wins, q_list, "白噪声 (期望单标度)")
    run_mfdfa(profile_mf, wins, q_list, "拼接序列 (近似多分形)")


if __name__ == "__main__":
    main()
​

💡 运行后把两条序列的 Δh、Δα 数字记下来，对照本章 5.6 节的判据表，亲手判定哪个是单标度、哪个是多分形。MFDFA 是 DFA 的自然推广：q=2 退化为 DFA；h(q) 谱宽揭示多分形性。

下一章分析两条序列的交叉标度：DCCA 与 MFDCCA。