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

本文介绍了一种稳定估计周期性平均跳跃点的方法，不依赖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便于验证(randomm在原代码中).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. 降噪:移动平均滤波(窗长需要窗长 > 噪声相关尺度，< K_true）
window = max(5, K_true // 3)  # 经验值：避免过度平滑，掩盖跳变
ma = np.convolve(y, np.ones(window)/window, mode='valid')
# 填充长度(前填充)
ma_padded = np.concatenate([np.full(window-1, ma[0]), ma])

# 3. 去趋势:拟合线性趋势并减去(因为均值阶梯整体呈负斜率)
z = np.polyfit(time, ma_padded, 1)
trend = np.poly1d(z)(time)
rm_trend = ma_padded - trend

# 4. 计算自相关函数（ACF）
corr = np.correlate(rm_trend - np.mean(rm_trend), 
                    rm_trend - np.mean(rm_trend), 
                    mode='full')
corr = corr[len(corr)//2:]  # 纠正滞后部分

# 5. 为ACF做FFT，识别主周期
freqs = fftfreq(len(corr), d=1.0)
fft_corr = fft(corr, norm='forward')
# 忽略零频(DC分量)，搜索非零频谱主峰
magnitudes = np.abs(fft_corr[1:len(corr)//2])
peak_idx, _ = find_peaks(magnitudes, height=np.max(magnitudes)*0.3)
if len(peak_idx) > 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()

关键注意事项：

选择窗口：移动平均窗口：移动平均窗口 window 应显着小于 K（建议 K//3 ~ K//2），否则会模糊阶跃边缘；
去除趋势的必要性：如果平均阶梯整体倾斜（如持续下降），必须去除线性趋势，否则 ACF 主峰会偏移；
K 整数约束：FFT 给出的是连续频率，需要四舍五入作为整数，并结合领域知识进行验证(如 K_est 应在 10–40 在合理范围内)；
替代方案：是的 rm_trend 直接使用 find_peaks(-rm_trend) 检测阶跃下降点，然后计算相邻峰间距，也可以稳定估计 K。

避免了这种方法 BIC 在信噪比中，模型假设负担充分利用了信号的周期性先验 ≥ 2 稳定的性能是工程实践中检测规则性平均跳变的首选策略。