变点检测 —— 一种贝叶斯方法
变点分析已经成为研究的许多领域的关注点。这种分析指的是在给定时间序列中找到突变或突然变化的问题。根据岩田等人（2018）的定义，变点分析是“识别时间序列发生概率分布变化的时刻的方法。”根据范登伯格和威廉斯（2020）的说法，“时间序列行为中突变的时刻通常是引起警报的原因，因为它们可能暗示数据生成过程发生了显著变化。”
正如阿米尼坎哈希和库克（2017）和岩田等人（2018）所示，对这种分析的增加关注是由于最近的技术发展。这些发展产生大量数据，通常需要密切监测，如机器人技术、医学、气象学、语音和图像识别等。有各种各样的模型和方法来处理这些问题。但由于本文不旨在对这些模型进行描述性分析，我想建议你查阅范登伯格和威廉斯（2020）的工作，以更好地了解这些方法。在那里，你可以找到在线和离线变点检测、单变量和多变量方法之间的区别，这些方法可能是参数化的或非参数化的，以及有监督或无监督的模型。
从一开始，根据埃勒斯（2007）的说法，时间序列“是随时间顺序观察的一系列观察结果”，其主要特征是给定观察对相邻观察的依赖性。时间序列可以是连续的或离散的，在前一种情况下，根据埃勒斯（2007）的说法，假设集合为T={t∶ t1< t < t2}，则时间序列表示为{X(t):t ∈ T}。因此，取时间序列T的观察窗口，其中有n个观察值，我们有一个由{Xm, X(m+1),…, Xn }表示的时间序列。根据阿米尼坎哈希和库克（2017，第3页）的说法，“变点的检测可以被定义为一个假设检验问题，涉及两个选择”，即“零假设H0：‘没有变化发生’‘和备择假设H1：‘发生变化’”。
开始
因此，如果你喜欢编码，现在是运行Jupyter Notebook并开始进行一些模拟和分析的时候了，因为我们将在这个方法中进行一次“随机漫步”。
让我们导入以下包：
import numpy as np
from numpy.random import seed
from numpy.random import randn
import random
import datetime
import matplotlib.pyplot as plt
import seaborn as sns
import math
import decimal
from scipy import stats
np.seterr(divide='raise') #Make sure you set this
我们将从一个包含两个时间序列的时间序列y(t)开始，其中y(t1)的均值为μ=1，加上一些噪声，以及y(t2)的均值为μ=2，加上一些噪声，两者都有30个观察值。正如你可能想象的那样，所提出的时间序列将有一个相当大的变点。
yt1 = u1 + 0.1*randn(30)
yt2 = u2 + 0.1*randn(30)
y = np.concatenate((yt1, yt2), axis=0)
plt.figure(figsize=(16, 10))
plt.plot(y)
如果你了解贝叶斯统计学，你会知道任何模型的构建基本上由3个分布组成。先验分布h(θ)反映了我们对问题的先前知识。似然函数f(x|θ)反映了获得的数据，并必须纳入先验分布。这将导致一个我们感兴趣的后验分布h(θ|x)。这就是我们使用贝叶斯定理的方式。
模糊聚类
在这一点上，我们面临的第一个（也是最重要的？）问题是从我们构建的时间序列中获得一个先验分布 —— 这是模型的第一部分。问题是：我们没有一个！如果你处理一个时间序列，一旦我们得到一个先验分布，大多数任务已经完成了。
达恩杰洛等人（2011）采用了一个有趣的方法来解决这个问题。他使用了Kohonen网络对时间序列进行聚类。与硬聚类不同，Kohonen网络是一种模糊聚类算法，这意味着任何给定点X都与概率p关联到A组。这种关联由函数fA(X)给出，该函数将A中的每个点与A中X的隶属度表示为区间[0, 1]内的实数。
有关Kohonen网络的完整且更好的解释，你可以参考Kohonen（1990）和Haykin（2007）。使用Python，我构建了这样一个网络，使用了两个函数：
def Center_Kohonen(y, X=0, K=2, alfa=0.8, C=500):
# This kohonen network implementation finds two centers in a time series
#Params:
#Y = Time Series
#M = Number of input patterns
#N = Dimension of the input pattern (for our case it will be 1)
#K = number of neurons, for the proposed problem (number of centers)
M = y.shape[0]
N = 1
f = 0
#Initializin the neurons' weights
I = y.argsort(axis=0) # Sorted Indexes
Y = np.sort(y) # Sortes Points in the Time Series
c1 = Y[0:7] # Beginning of the series
c2 = Y[M-7:M] # End of the series
#Adjusting the values
while np.std(c2) > 0.1: #As long as the standard deviation is greater than 0.1, replace the highest value with the mean
ma = c2.argmax()
c2[ma] = np.mean(c2);
y[I[ma+60-7]] = np.mean(c2);
while np.std(c1) > 0.1: #As long as the standard deviation is greater than 0.1, replace the lowest value with the mean
mi = c1.argmin()
c1[mi] = np.mean(c1)
y[I[mi]] = np.mean(c1);
#Definition of weight values
W = [np.mean(c1), np.mean(c2)]
#Finding centers from Kohonen's network training
for l in range(1, C+1):
alfa=alfa*(1-(C-(C-l))/C)
for i in range(0, M): #For each value in the time series
a=999
for j in range(0, K): #Where K is the number of cluster
if np.linalg.norm([(y[i]-W[j])], 2) < a:
a = np.linalg.norm([(y[i]-W[j])], 2)
f = j
W[f] = W[f]+alfa*(y[i]-W[f])
return c1, c2, I, Y, W, a, alfa
def Fuzzy_Set(y, W):
# This program finds membership values for a time series from previously specified centers.
# Where y is the time series and c a vector with the found centers
center_1 = []
center_2 = []
n = y.shape[0]
l = 2
# Finding the membership association for each point in the time series
for i in range(0, l):
for t in range(0, n):
sum=0;
for k in range(0, l):
sum = sum+(y[t]-W[k])*(y[t]-W[k])
if i == 0:
center_1.append(np.round(1-((y[t]-W[i])*(y[t]-W[i]))/sum, 3))
else:
center_2.append(np.round(1-((y[t]-W[i])*(y[t]-W[i]))/sum, 3))
return center_1, center_2
如果你依次调用这两个函数并绘制结果，你可能会得到类似于这样的图表：
c1, c2, I, Y, W, a, alfa = Center_Kohonen(y, X=0, K=2, alfa=0.8, C=500)
center_1, center_2 = Fuzzy_Set(y, W)
plt.figure(figsize=(16, 10))
plt.plot(center_1, 'b') # plotting t, a separately
plt.plot(center_2, 'r') # plotting t, b separately
plt.show()
这非常有趣！通过Kohonen网络，我们能够将时间序列y(t)拆分。这张图表显示我们有两个聚类，因为我们设置了K=2。在横轴上，我们有时间序列的每个点，而纵轴显示了给定点与两个聚类之一关联的概率。正如你可能看到的，蓝线显示我们在达到第30点之前，所有点更有可能（约99%左右）属于第一组或集合μ1(t)。红线显示相反的情况，因为它表示与第二组即集合μ2(t)的关联。这是合理的，因为我们构建了一个具有两个不同均值的时间序列，而且形象地说，这个图与第一个图相关。
尽管有趣，但直到现在我们并没有真正找到变点（我们有一些线索），而且这里没有贝叶斯的内容。
顺便说一下，大多数时候，时间序列的点之间的区分不会那么简单。例如，如果我们构建了时间序列y(t)，其中y(t1)的均值为μ=1加上一些噪声，而y(t2)的均值为μ=1.3（而不是2）加上一些噪声，那么拆分会好吗？我让你尝试一下这个练习……
Metropolis-Hastings算法登场
如果你尝试了上面的练习，你会发现只使用Kohonen网络来找到替代时间序列中任何变点的迹象会让你陷入困境。这是因为Kohonen网络并不提供变点，而是提供两组连续变量，表示每个点与给定聚类的关联。
但请记住，集合μ1(t)和μ2(t)的值在区间[0,1]内。这意味着μ1(t)和μ2(t)是具有不同参数的Beta分布的近似值（你听说过Kullback–Leibler吗？）。根据达恩杰洛等人（2011）的说法，假设变点由m表示，那么对于t≤m，我们将有一个Beta(a,b)分布，对于t>m，我们将有一个Beta(c,d)分布。根据Beta分布的性质，如果时间序列中有任何变点，在Beta(a,b)中参数a将大于参数b，而在Beta(c,d)中参数c将小于参数d。
问题是：你如何构建这两个Beta分布？Metropolis-Hastings算法是一种马尔可夫链蒙特卡洛的一种形式，最初由Metropolis等人（1953）提出，然后由Hastings（1970）进行了推广。根据Gelman等人（2003）的说法，任何马尔可夫链模拟的目标“是创建一个稳态分布由p(θ | x)指定的马尔可夫过程。”运行足够的模拟可以使我们获得足够接近稳态和后验分布的分布。后验分布可以用参数θ的某些函数的期望来总结，即∫g(θ)p(θ | x)dθ = E [g(θ) | x]。这种积分并不是微不足道的，这就是MCMC方法用于近似后验分布的原因。
Metropolis-Hastings算法使用拒绝的思想，这意味着它从辅助分布生成一个值，并以给定的概率接受它。如果你对MCMC方法不熟悉，你可能会质疑算法如何拒绝抽取的值。我们使用Hastings（1970）给出的转移规则：
换句话说，我们可以使用由模糊聚类给出的两组连续变量，从给定变点检测的先验分布中随机拒绝抽取的值。如果你想了解更多关于MCMC方法的信息，我建议参考Gamerman和Lopes（2018）。
让我们回到Jupyter Notebook。下面的函数是为这个问题实现Metropolis-Hastings算法。尽管功能强大，但该算法需要一些东西。首先是为需要找到的每个参数设置一个先验分布。对于参数m，我们使用1到60之间的均匀分布，这意味着算法在时间序列中随机选择一个变点候选。对于参数a、b、c和d，我选择了弱信息的伽马分布。该函数还需要参数，这是一组随机变量（μ1(t)或μ2(t)）和模拟次数。
def Metropolis_Hastings(center_kohonen, n_sims=1000):
n = len(y)
m = 1 + round((n-1) * np.random.uniform(0, 1))
shape, scale, loc = 10, 0.1, 0
#Lists to save the date for each parameter
a_params = []
b_params = []
c_params = []
d_params = []
m_params = []
#Prior Distributions for the Parameters
a = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
b = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
c = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
d = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
for i in range(0, n_sims):
m1 = 1+round((n-1) * np.random.uniform(0, 1));
a1 = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
b1 = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
c1 = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
d1 = stats.erlang.rvs(shape, loc=loc, scale=scale, size=1, random_state=None)[0]
#PARAM A
aux1 = 1
for j in range(0, m):
try:
aux1 = aux1 * (center_kohonen[j] ** (a1-1))
except:
aux1 = aux1
aux2 = 1
for j in range(0, m):
try:
aux2 = aux2 * center_kohonen[j] ** (a-1)
except:
aux2 = aux2
try:
ra = ((math.gamma(a1+b)/math.gamma(a1))**m)*aux1*((((a/a1)**.9)*math.exp(-.1*(a1-a)))**2)/(((math.gamma(a+b)/math.gamma(a))**m)*aux2)
if (min(1, ra) > np.random.uniform(0, 1)):
a=a1
except:
pass
#PARAM B
aux1 = 1
for j in range(0, m):
try:
aux1 = aux1*(1-center_kohonen[j])**(b1-1)
except:
aux1 = aux1
aux2 = 1
for j in range(0, m):
try:
aux2 = aux2*(1-center_kohonen[j])**(b-1)
except:
aux2 = aux2
try:
rb = ((math.gamma(a+b1)/math.gamma(b1))**m)*aux1*((((b/b1)**.9)*math.exp(-.1*(b1-b)))**2)/(((math.gamma(a+b)/math.gamma(b))**m)*aux2)
if (min(1, rb) > np.random.uniform(0, 1)):
b = b1
except:
pass
#PARAM C
aux1 = 1
for j in range(m, n):
try:
aux1=aux1*center_kohonen[j]**(c1-1)
except:
aux1 = aux1
aux2 = 1
for j in range(m, n):
try:
aux2=aux2*center_kohonen[j]**(c-1)
except:
aux2 = aux2
try:
rc = ((math.gamma(c1+d)/math.gamma(c1))**(n-m))*aux1*((((c/c1)**.9)*math.exp(-.1*(c1-c)))**2)/(((math.gamma(c+d)/math.gamma(c))**(n-m))*aux2)
if (min(1, rc) > np.random.uniform(0, 1)):
c = c1
except:
pass
#PARAM D
aux1 = 1
for j in range(m, n):
try:
aux1=aux1*(1-center_kohonen[j])**(d1-1)
except:
aux1 = aux1
aux2 = 1
for j in range(m, n):
try:
aux2=aux2*(1-center_kohonen[j])**(d-1)
except:
aux2 = aux2
try:
rd = ((math.gamma(c+d1)/math.gamma(d1))**(n-m))*aux1*((((d/d1)**.9)*math.exp(-.1*(d1-d)))**2)/(((math.gamma(c+d)/math.gamma(d))**(n-m))*aux2)
if (min(1, rd) > np.random.uniform(0, 1)):
d = d1
except:
pass
#PARAM M
aux1 = 1
for j in range(0, m1):
try:
aux1 = aux1*(center_kohonen[j]**(a-1))*((1-center_kohonen[j])**(b-1))
except:
aux1 = aux1
aux2 = 1;
for j in range(m1, n):
try:
aux2 = aux2*(center_kohonen[j]**(c-1))*((1-center_kohonen[j])**(d-1))
except:
aux2 = aux2
aux3 = 1
for j in range(0, m):
try:
aux3 = aux3*(center_kohonen[j]**(a-1))*((1-center_kohonen[j])**(b-1))
except:
aux3 = aux3
aux4 = 1
for j in range(m, n):
try:
aux4 = aux4*(center_kohonen[j]**(c-1))*((1-center_kohonen[j])**(d-1))
except:
aux4 = aux4
try:
rm = (((math.gamma(a+b)/(math.gamma(a)*math.gamma(b)))**m1)*((math.gamma(c+d)/(math.gamma(c)*math.gamma(d)))**(n-m1))*aux1*aux2)/(((math.gamma(a+b)/(math.gamma(a)*math.gamma(b)))**m)*((math.gamma(c+d)/(math.gamma(c)*math.gamma(d)))**(n-m))*aux3*aux4)
if (min(1, rm) > np.random.uniform(0, 1)):
m = m1
except:
pass
a_params.append(a)
b_params.append(b)
c_params.append(c)
d_params.append(d)
m_params.append(m)
return a_params, b_params, c_params, d_params, m_params
使用两个必需参数调用该函数：在这里，我发送了由函数Fuzzy_Set给出的center_1和n_sims=1000
a_params, b_params, c_params, d_params, m_params = Metropolis_Hastings(center_1, n_sims=1000)
fig_dims = (16, 10)
fig, ax = plt.subplots(figsize=fig_dims)
plt.plot(m_params, 'r')
ax.set(xlabel='# Simulations', ylabel='Change Point Candidates (m)')
现在你终于找到了变点。这张图表很有趣，因为它显示了抽取过程是如何进行的。由均匀分布给出的第一个值是m=55。算法拒绝了它，然后尝试另一个，直到获得令人满意且稳定的结果。经过约150次额外运行后，m=30的值不再被算法拒绝。 由于该函数返回了每个参数的采样值，我们也可以绘制它们的值。从参数m开始，它是所有抽取的变点。为了查看密度图，你可以舍弃前200次模拟作为“burn-in”：
fig_dims = (16, 10)
fig, ax = plt.subplots(figsize=fig_dims)
ax.set(xlabel='Change Point Candidates', ylabel='Density')
sns.kdeplot(m_params[200:])
我们还可以使用参数a、b、c和d的平均值创建Beta分布。正如我们之前讨论的，这些参数对于Metropolis-Hastings算法至关重要，因为拒绝规则必须断言对于t≤m，我们将有一个Beta(a,b)分布，对于t>m，我们将有一个Beta(c,d)分布。
让我们继续构建这样的表示，使用包含参数a、b、c和d的采样值的变量a_params、b_params、c_params和d_params：
fig, ax = plt.subplots(1, 1, figsize=(16, 8))
ax.set_ylabel('Density')
beta_1 = np.random.beta(np.mean(a_params[200:]), np.mean(b_params[200:]), size=1000)
beta_2 = np.random.beta(np.mean(c_params[200:]), np.mean(d_params[200:]), size=1000)
sns.kdeplot(beta_1, color='b', fill=True, common_norm=False, alpha=.5, linewidth=0)
sns.kdeplot(beta_2, color='r', fill=True, common_norm=False, alpha=.5, linewidth=0)
用参数a和b的平均值构建的第一个Beta分布是红色的，而用参数c和d的平均值构建的第二个是蓝色的。在中间，两者都不太密集的地方，我们找到了变点。无疑，使用这种方法的一个巨大优势是可以获得这样的分布，因为我们可以用它们进行贝叶斯推断，丰富预测模型，甚至使用其他类型的蒙特卡洛模拟。
结论
在时间序列中找到变点可以防止系统陷入严重的麻烦。想象一下一个需要控制温度的机器。任何突然的变化都必须尽快被识别，以便工程师检查它们。或者在各种生产设施中的能源消耗。任何过度消耗都必须被分析，因为它可能表示生产中的某些偏差，甚至能源泄漏，从而显著影响生产成本。
话虽如此，D'Angelo（2011）在这里在Python中实现的方法被证明对检测给定时间序列中的变点非常有价值。此外，正如先前指出的，这种方法的另一个巨大优势正是我们获得了两个Beta分布作为输出，这可能非常有用。
参考资料
AMINIKHANGHAHI, Samaneh e COOK, J. Diane. A Survey of Methods for Time Series Change Point Detection. Knowledge and Information Systems, 51, 2017.
D’ANGELO, Marcos Flávio et al. A Fuzzy/Bayesian Approach for the Time Series Change Point Detection Problem. Pesquisa Operacional, 31(2), 2011.
EHLERS, Ricardo S. Inferência Bayesiana. 5ª Ed. Departamento de Estatística da Universidade Federal do Paraná, 2007.
FAMA, Eugene. Random Walks in Stock Market Prices. Financial Analysts, 21, 1965.
_ Efficient Capital Markets: A Review of Theory and Empirical Work. The Journal of Finance, 25, 1970.
GAMERMAN, Dani. LOPES, Hedibert Freitas. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. 2º Ed. Florida : Chapman & Hall/CRC, 2006. 315p.
GELMAN, Andrew et al. Bayesian Data Analysis. 2ª Ed. Florida : Chapman & Hall/CRC, 2003. 668p.
HASTINGS, W. K. Monte Carlo Sampling Methods Using Markov Chains and their Applications. Biometrika, 57, 1970.
HAYKIN, Simon. Redes Neurais: Princípios e Práticas. 2ª Ed. Porto Alegre : Bookman, 2007. 900p.
IWATA, Takuma et al. Accelerating Online Change-Point Detection Algorithm using 10GbE FPGA NIC. Trabalho apresentado na 24ª Conferência Internacional de Computação Paralela e Distribuída. Turim : 2018.
KOHONEN, Teuvo. Self-Organizing Maps. Nova Iorque : Springer, 1990, 501p.
METROPOLIS, Nicholas et al. Equation of State Calculations by Fast Computing Machines. Journal of Chemical Physics, 21, 1953.
OH KJ, et al. Developing Time-Based Clustering Neural Networks To Use Change-Point Detection: Application To Financial Time Series. Asia-Pacific Journal of Operational Research, 22(1), 2005.
PAULINO, Carlos Daniel et al. Estatística Bayesiana. 2º Edição. Lisboa : Fundação Calouste Gulbenkian, 2018, 601p.
VAN DEN BURG, Gerrit J. J., WILLIAMS, Christopher K. I. An Evaluation of Change Point Detection Algorithms. stat.ML, arXiv:2003.06222, 2020
✄-----------------------------------------------