示例11 分式布朗运动合成与标定

"""ex11_fbm_synthesis.py — 第 10 章 合成分式布朗运动 (fBm) 并标定 Hurst

目标：用 Davies-Harte 近似方法生成给定 Hurst 的分式布朗运动 (fBm)，
然后用 DFA 反推 H，验证 fathon 在已知答案的合成数据上的准确度。

对应教程：tutorials/10-advanced.md。

运行：
python ex11_fbm_synthesis.py
"""

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

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

fBm 的自协方差函数：
    γ(k) = 0.5 * (|k+1|^(2H) - 2|k|^(2H) + |k-1|^(2H))

利用循环嵌入 + FFT，得到长度为 n 的 fBm 样本。
fBm 本身已经是 profile-like（累积和）序列。
"""
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))
# 循环嵌入：长度 2n，第二半镜像
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]
# Davies-Harte 的输出量纲需要归一化到标准差
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}")
results = []
for H_true in [0.3, 0.5, 0.7, 0.9]:
    # fBm 的「增量」才是 fGn；DFA 的 toAggregated 等价于再累积一次
    # 教程惯例：直接对 fBm 样本做 toAggregated（或对 fGn 做）
    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}")
    results.append((H_true, H_est))

# ── 可视化不同 H 的 fBm 形态 ──
fig, axes = plt.subplots(2, 2, figsize=(11, 6), sharex=True)
for ax, H_true in zip(axes.flat, [0.3, 0.5, 0.7, 0.9]):
    bm = fbm_davies_harte(H_true, 1000, seed=int(H_true * 100))
    ax.plot(bm, lw=0.8)
    ax.set_title(f'H = {H_true:.1f}  '
                 f'{"反持久" if H_true < 0.5 else "白噪声" if H_true == 0.5 else "持久"}')
    ax.grid(alpha=0.3)
plt.suptitle('不同 Hurst 指数的分式布朗运动形态', fontsize=13)
plt.tight_layout()
plt.show()
​

if name == "main":
main()