时间序列分散的步数近似，使用 BIC 改变每 K 个步数的均值

Question

首先，所使用的合成数据是按以下方式生成的：

import sympy as sp 
import numpy as np
import matplotlib.pyplot as plt
import random
import math

np.random.seed(2)

n_samples = 180
time = np.arange(n_samples)

mean_value = random.randrange(60, 90)
mean = np.full(n_samples, mean_value)

# The fixed mean segment is generated randomly
K = random.randint(10, 40)

for i in range(K, n_samples, K):
    mean[i:] = mean[i - K] - 10  

noise = np.random.randn(n_samples) * random.normalvariate(4, 2)
y = mean + noise

我不知道如何近似均值并检测均值的变化，考虑到涉及噪声并且方差在各个步骤中是恒定的，但仍然未知，到目前为止，我将似然函数 L 作为正态似然函数，但我不知道如何使用它知道 BIC = -2Log(L) 到目前为止我的代码是

def find_optimal_change_point(data):
    min_bic = float('inf')
    bics = np.full(len(data),min_bic)
    change_points = [0] * K

    for i in range(1, len(data)):  
        segment = data[:i]
        
        mean, var = np.mean(data[:i - 1]), np.var(data[:i - 1])

        N = len(segment)
        S = var
        new_bic = bic(N, S, 4, segment, v, len(segment),index=i) 

        if bics[i-1] > new_bic:
            bics[i] = new_bic
            change_points[k] = data[i]
            
        elif bics[i] < bics[i-1] :
            break
        
    return change_points

我需要在散点图上绘制步骤，其中步骤中的垂直线遵循每个步骤的平均值，步骤的高度连接平均值发生变化的两个点

Answer 1

我不确定您为什么要使用 BIC，因为它通常是模型选择的工具，而不是有意义的时间序列统计数据。

另一种效果很好的方法可能是平滑噪声信号以消除噪声（例如通过移动平均值），并消除信号中的趋势。然后使用傅里叶变换和/或相关性来检测信号的周期性（这应该是均值变化的周期）。 由此应该很容易得出平均值。

这是我测试过的一个小例子，作为第一个近似值效果很好：

from scipy.fft import fft, fftfreq
n = 50
ma = np.convolve(y,np.ones(n), mode='valid')/n # denoised signal
rm_trend = y-((ma[-1]-ma[0])/len(ma)*np.arange(len(y))+ma[0]) # remove trend

corr = np.correlate(rm_trend,rm_trend,mode='full')
corr = corr[corr.shape[0]//2:]

y_fft = fft(y,norm='forward')[1:len(y)//2] # remove the mean

corr = np.correlate(rm_trend,rm_trend,mode='full') # autocorrelation
corr = corr[corr.shape[0]//2:]
freq = fftfreq(len(corr)) # frequencies
corr_fft = fft(corr,norm='forward')[1:len(corr)//2] # FFT without mean
k = 1/freq[np.argmax(corr_fft )+1]
print(k)

如果这不能回答您的问题，请告诉我。