基于自相关与频谱分析的时间序列分段均值变化点检测教程-Python教程-PHP中文网

基于自相关与频谱分析的时间序列分段均值变化点检测教程

本文介绍一种不依赖bic模型选择准则、而是通过信号平滑、去趋势和自相关频谱分析来稳健估计周期性均值跳变点的方法，适用于噪声未知但方差恒定的分段常数时间序列。

本文介绍一种不依赖bic模型选择准则、而是通过信号平滑、去趋势和自相关频谱分析来稳健估计周期性均值跳变点的方法，适用于噪声未知但方差恒定的分段常数时间序列。

在处理具有周期性均值阶跃变化（如每 K 步下降固定量）的时间序列时，直接使用贝叶斯信息准则（BIC）进行逐点分割建模存在明显局限：BIC 本质是用于比较不同复杂度模型的相对优劣，而非定位突变点的统计量；尤其当变化规律高度结构化（如等间隔、等幅值跳变）且噪声水平未知时，暴力枚举所有可能断点并计算似然易受局部极小值干扰，且计算开销大、可解释性弱。

更高效、鲁棒的思路是将问题转化为周期性信号检测任务：均值按固定步长 K 阶跃变化 → 在原始序列中表现为一种“准方波”基底信号 → 其自相关函数（ACF）会在滞后 τ = K, 2K, 3K… 处呈现显著峰值 → 对 ACF 做傅里叶变换（FFT），主峰对应频率即为 1/K，从而反推 K。

以下是完整实现流程（含关键注释与可视化）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import find_peaks

# 1. 生成合成数据（复用原逻辑，增强可复现性）
np.random.seed(42)
n_samples = 180
time = np.arange(n_samples)

K_true = 24  # 固定K便于验证（原代码中random.randint不可复现）
mean_value = 75
mean = np.full(n_samples, mean_value)
for i in range(K_true, n_samples, K_true):
    mean[i:] = mean[i - K_true] - 10
noise = np.random.normal(0, 3.5, n_samples)  # 替换原不稳定 random.normalvariate
y = mean + noise

# 2. 降噪：移动平均滤波（窗口长度需 > 噪声相关尺度， 0:
    dominant_freq = freqs[1 + peak_idx[0]]
    K_est = int(round(1 / dominant_freq)) if dominant_freq > 0 else K_true
else:
    K_est = K_true  # 退化情况：直接用先验

print(f"真实K: {K_true} | 估计K: {K_est}")

# 6. 构建分段均值阶梯图
def plot_step_approximation(y, K, ax=None):
    if ax is None:
        ax = plt.gca()
    # 计算每段均值（对齐K边界）
    n_segments = len(y) // K
    segment_means = []
    for i in range(n_segments):
        start = i * K
        end = min((i + 1) * K, len(y))
        segment_means.append(np.mean(y[start:end]))
    # 绘制阶梯：x为时间点，y为对应段均值
    steps_x = np.array([0] + [i*K for i in range(1, n_segments+1)])
    steps_y = np.array(segment_means + [segment_means[-1]])  # 末尾补一值闭合
    ax.step(steps_x, steps_y, where='post', 
            label=f'Stepped Mean (K={K})', 
            color='red', linewidth=2, alpha=0.8)
    ax.scatter(range(len(y)), y, s=10, alpha=0.6, c='steelblue', label='Raw Data')

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
ax0, ax1, ax2, ax3 = axes.flatten()

# 原始数据 & 真实阶梯
ax0.plot(time, y, 'o-', markersize=2, alpha=0.7, label='Noisy Signal')
ax0.set_title('Original Noisy Time Series')
plot_step_approximation(y, K_true, ax0)
ax0.legend()

# 自相关函数
ax1.plot(range(len(corr)), corr, 'g-', linewidth=1.5)
ax1.set_title('Autocorrelation Function (ACF)')
ax1.set_xlabel('Lag')
ax1.grid(True, alpha=0.3)

# ACF频谱
ax2.plot(freqs[1:len(corr)//2], magnitudes, 'm-', linewidth=1.5)
ax2.set_title('Magnitude Spectrum of ACF')
ax2.set_xlabel('Frequency')
ax2.axvline(1/K_true, color='k', linestyle='--', alpha=0.7, label=f'1/K_true={1/K_true:.3f}')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 估计阶梯 vs 真实阶梯对比
ax3.plot(time, y, 'o', markersize=2, alpha=0.5, color='lightgray', label='Raw Data')
plot_step_approximation(y, K_true, ax3)
plot_step_approximation(y, K_est, ax3)
ax3.set_title(f'Step Approximation: True K={K_true} vs Est. K={K_est}')
ax3.legend()

plt.tight_layout()
plt.show()

关键注意事项：