示例05 MFDFA 多分形四联图

"""ex05_mfdfa.py — 第 5 章 MFDFA 多分形扩展

目标：跑通 MFDFA，得到广义 Hurst h(q)、质量指数 τ(q)、多重分形谱 f(α)，
并用四联图直观展示「单标度」与「多分形」的差别。

对应教程：tutorials/05-mfdfa.md。

运行：
python ex05_mfdfa.py
"""

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

def make_multifractal_like(n, seed=0):
"""合成一个近似多分形的序列：不同尺度分段白噪声拼接，
使 h(q) 随 q 变化（教程 5.4 节思路）。"""
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()
print(f"\n[{title}]")
print(f" h(2) = {h_q[np.argmin(np.abs(q_list - 2))]:.4f} (应对标 DFA 的 H)")
print(f" Δh = {h_q.max() - h_q.min():.4f}")
print(f" Δα = {alpha.max() - alpha.min():.4f}")
return dict(n=n, F=F, h_q=h_q, tau=tau, alpha=alpha, mfSpect=mfSpect,
q_list=q_list, title=title)

def plot_quad(res):
n, F, h_q, tau, alpha, mfSpect, q_list, title = (
res['n'], res['F'], res['h_q'], res['tau'],
res['alpha'], res['mfSpect'], res['q_list'], res['title']
)
fig, axes = plt.subplots(2, 2, figsize=(11, 9))

# 1. F_q(n) 标度
for qi in [0, 2, 4, -2]:
    idx = int(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})')

fig.suptitle(title, fontsize=13)
plt.tight_layout()
​

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)

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

plot_quad(res_white)
plot_quad(res_mf)
plt.show()

print("\n判据：Δh < 0.1 视为单标度；> 0.1 视为多分形。")
​

if name == "main":
main()