R中时间序列中的变化点检测

Question

我需要一些关于变更点如何在时间序列中工作的指导。我正在尝试使用 R 检测一些变更点，以及名为“changepoint”的包 (https://cran.r-project.org/web/packages/changepoint/changepoint.pdf)。

有一些选项可用于检测方差 (cpt.var) 和均值 (cpt.mean) 何时发生变化，但我要寻找的是时间序列何时变化趋势。

也许我对真正的变更点感到困惑，但有什么方法可以获取这些信息？

我正在显示使用 cpt.var() 函数的结果，并且我添加了一些箭头，显示了我想要实现的目标。

有什么方法可以实现吗？我想应该有点像拐点......

我将不胜感激。

先谢谢了， 乔恩

编辑

我尝试过使用diff() 的方法，但没有正确检测到更改：

我使用的数据如下：

  [1] 10.695 10.715 10.700 10.665 10.830 10.830 10.800 11.070 11.145 11.270 11.015 11.060 10.945 10.965 10.780 10.735 10.705 10.680 10.600 10.335 10.220 10.125
 [23] 10.370 10.595 10.680 11.000 10.980 11.065 11.060 11.355 11.445 11.415 11.350 11.310 11.330 11.360 11.445 11.335 11.275 11.300 11.295 11.470 11.445 11.325
 [45] 11.300 11.260 11.200 11.210 11.230 11.240 11.300 11.250 11.285 11.215 11.260 11.395 11.410 11.235 11.320 11.475 11.470 11.685 11.740 11.740 11.700 11.905
 [67] 11.720 12.230 12.285 12.505 12.410 11.995 12.110 12.005 11.915 11.890 11.820 11.730 11.700 11.660 11.685 11.615 11.360 11.425 11.185 11.275 11.265 11.375
 [89] 11.310 11.250 11.050 10.880 10.775 10.775 10.805 10.755 10.595 10.700 10.585 10.510 10.290 10.255 10.395 10.290 10.425 10.405 10.365 10.010 10.305 10.185
[111] 10.400 10.700 10.725 10.875 10.750 10.760 10.905 10.680 10.670 10.895 10.790 10.990 10.925 10.980 10.975 11.035 10.895 10.985 11.035 11.295 11.245 11.535
[133] 11.510 11.430 11.450 11.390 11.520 11.585

当我执行 diff() 时，我会得到以下数据：

  [1]  0.020 -0.015 -0.035  0.165  0.000 -0.030  0.270  0.075  0.125 -0.255  0.045 -0.115  0.020 -0.185 -0.045 -0.030 -0.025 -0.080 -0.265 -0.115 -0.095  0.245
 [23]  0.225  0.085  0.320 -0.020  0.085 -0.005  0.295  0.090 -0.030 -0.065 -0.040  0.020  0.030  0.085 -0.110 -0.060  0.025 -0.005  0.175 -0.025 -0.120 -0.025
 [45] -0.040 -0.060  0.010  0.020  0.010  0.060 -0.050  0.035 -0.070  0.045  0.135  0.015 -0.175  0.085  0.155 -0.005  0.215  0.055  0.000 -0.040  0.205 -0.185
 [67]  0.510  0.055  0.220 -0.095 -0.415  0.115 -0.105 -0.090 -0.025 -0.070 -0.090 -0.030 -0.040  0.025 -0.070 -0.255  0.065 -0.240  0.090 -0.010  0.110 -0.065
 [89] -0.060 -0.200 -0.170 -0.105  0.000  0.030 -0.050 -0.160  0.105 -0.115 -0.075 -0.220 -0.035  0.140 -0.105  0.135 -0.020 -0.040 -0.355  0.295 -0.120  0.215
[111]  0.300  0.025  0.150 -0.125  0.010  0.145 -0.225 -0.010  0.225 -0.105  0.200 -0.065  0.055 -0.005  0.060 -0.140  0.090  0.050  0.260 -0.050  0.290 -0.025
[133] -0.080  0.020 -0.060  0.130  0.065

我得到的是下一个结果：

> cpt =cpt.mean(diff(vector), method="PELT")

> (cpt.pts <- attributes(cpt)$cpts)
[1] 137

似乎没有意义...有什么线索吗？

Answer 1

如果信号不是太嘈杂，您可以使用diff 来检测斜率变化点而不是平均值：

library(changepoint)

set.seed(1)
slope <- rep(sample(10,10)-5,sample(100,10))
sig <- cumsum(slope)+runif(n=length(slope),min = -1, max = 1)
cpt =cpt.mean(diff(sig),method="PELT")

# Show change point
(cpt.pts <- attributes(cpt)$cpts)
#> [1]  58 109 206 312 367 440 447 520 599

plot(sig,type="l")
lines(x=cpt.pts,y=sig[cpt.pts],type="p",col="red",cex=2)

另一个似乎更适合您提供的数据的选项是使用piecewise linear segmentation：

library(ifultools)
changepoints <- linearSegmentation(x=1:length(data),y=data,angle.tolerance = 90,n.fit=10,plot=T)
changepoints
#[1]  13  24  36  58  72 106

Answer 2

在 R 中，有许多包可用于时间序列变化点检测。 changepoint 绝对是一个非常有用的。 CRAN 任务视图中总结了部分软件包列表：

结构变化（使用线性回归模型）和趋势（使用非参数测试）中提供了变化点检测。 changepoint 包提供了许多流行的 changepoint 方法，ecp 对单变量和多变量序列进行非参数化的 changepoint 检测。 changepoint.np 实现了非参数 PELT 算法，而 changepoint.mv 检测多元时间序列中的变化点。 InspectChangepoint 使用稀疏投影来估计高维时间序列中的变化点。 robcp 使用 Huberized cusum 测试提供稳健的变化点检测，而 Rbeast 提供贝叶斯变化点检测和时间序列分解。

这里还有一个很棒的博客，比较了几个替代包：https://www.marinedatascience.co/blog/2019/09/28/comparison-of-change-point-detection-methods/。另一个令人印象深刻的比较来自 Jonas Kristoffer Lindeløv 博士，他开发了 mcp 软件包：https://lindeloev.github.io/mcp/articles/packages.html。

在下面，我使用您的示例时间序列使用我自己开发的 Rbeast 包生成了一些快速结果（选择这里显然是为了自我推销和感知相关性）。 Rbeast是一种贝叶斯变化点检测算法，可以估计变化点发生的概率。它还可以用于将时间序列分解为季节性和趋势，但显然，您的时间序列是仅趋势的，因此在下面的beast 函数中，指定了season='none'。

y = c(10.695,10.715,10.700,10.665,10.830,10.830,10.800,11.070,11.145,11.270,11.015,11.060,10.945,10.965,10.780,10.735,10.705,
    10.680,10.600,10.335,10.220,10.125,10.370,10.595,10.680,11.000,10.980,11.065,11.060,11.355,11.445,11.415,11.350,11.310,11.330,
    11.360,11.445,11.335,11.275,11.300,11.295,11.470,11.445,11.325,11.300,11.260,11.200,11.210,11.230,11.240,11.300,11.250,11.285,
    11.215,11.260,11.395,11.410,11.235,11.320,11.475,11.470,11.685,11.740,11.740,11.700,11.905,11.720,12.230,12.285,12.505,12.410,
    11.995,12.110,12.005,11.915,11.890,11.820,11.730,11.700,11.660,11.685,11.615,11.360,11.425,11.185,11.275,11.265,11.375,11.310,
    11.250,11.050,10.880,10.775,10.775,10.805,10.755,10.595,10.700,10.585,10.510,10.290,10.255,10.395,10.290,10.425,10.405,10.365,
    10.010,10.305,10.185,10.400,10.700,10.725,10.875,10.750,10.760,10.905,10.680,10.670,10.895,10.790,10.990,10.925,10.980,10.975,
    11.035,10.895,10.985,11.035,11.295,11.245,11.535 ,11.510,11.430,11.450,11.390,11.520,11.585)

library(Rbeast)
out=beast(y, season='none')
plot(out)
print(out)

在上图中，垂直虚线标记了最可能的变化点位置； Pr(tcp) 的绿色曲线显示随时间发生变化点的逐点概率。 order_t 曲线给出了充分拟合趋势所需的分段多项式的估计平均阶（第 0 阶是恒定的，第 1 阶是线性的）：趋于 0 的平均阶意味着趋势更可能是平坦的，并且接近 1 的订单意味着趋势是线性的。输出也可以打印为一些 ascii 输出，如下所示。同样，它说时间序列最有可能有 8 个变化点；他们最可能的位置在out$trend$cp 中给出。

Result for time series #1 (total number of time series in 'out': 1)

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+                     SEASONAL CHANGEPOINTS                    +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


 No seasonal/periodic component present (i.e., season='none')


++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+                     TREND CHANGEPOINTS                       +
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


An ascii plot of the probability dist for number of chgpts(ncp)
---------------------------------------------------------------
Pr(ncp=0 )=0.000|*                                            |
Pr(ncp=1 )=0.000|*                                            |
Pr(ncp=2 )=0.000|*                                            |
Pr(ncp=3 )=0.000|*                                            |
Pr(ncp=4 )=0.000|*                                            |
Pr(ncp=5 )=0.000|*                                            |
Pr(ncp=6 )=0.055|*****                                        |
Pr(ncp=7 )=0.074|******                                       |
Pr(ncp=8 )=0.575|******************************************** |
Pr(ncp=9 )=0.240|*******************                          |
Pr(ncp=10)=0.056|*****                                        |
---------------------------------------------------------------
Max ncp : 10   | A parameter you set (e.g., maxTrendKnotNum)  |
Mode ncp: 8    | Pr(ncp= 8)=0.57; there is a 57.5% probability|
           | that the trend componet has  8 chngept(s).   |
Avg ncp : 8.17 | Sum[ncp*Pr(ncp)]                             |
---------------------------------------------------------------

List of most probable trend changepoints (avg number of changpts: 8.17) 
--------------------------------.
tcp# |time (cp)      |prob(cpPr)|
-----|---------------|----------|
1    |8.0000         |   0.92767|
2    |112.0000       |   0.91433|
3    |68.0000        |   0.84213|
4    |21.0000        |   0.80188|
5    |32.0000        |   0.78171|
6    |130.0000       |   0.76938|
7    |101.0000       |   0.66404|
8    |62.0000        |   0.61171|
--------------------------------'