R中时间序列中的变化点检测
我需要一些关于转换点如何在时间序列中工作的指导。我正在尝试使用R来检测一些转换点，以及称为“转换点”(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阿佩莉没有道理..。有线索吗？
回答 2
Stack Overflow用户
发布于 2022-01-15 05:30:36
在R中，有许多软件包可用于时间序列转换点的检测。changepoint绝对是非常有用的。在CRAN任务视图中总结了包的部分列表：
变化点检测是在结构变化(使用线性回归模型)和趋势(使用非参数检验)中提供的。转换点包提供了许多流行的转换点方法，ecp为单变量和多元序列进行非参数转换点检测。changepoint.np实现了非参数PELT算法，而changepoint.mv在多变量时间序列中检测转换点。InspectChangepoint使用稀疏投影来估计高维时间序列中的转换点.robcp利用Huberized测试提供了鲁棒的变化点检测，Rbeast提供了贝叶斯变化点检测和时间序列分解。
这里还有一个很棒的博客，比较了几个可供选择的软件包：https://www.marinedatascience.co/blog/2019/09/28/comparison-of-change-point-detection-methods/。另一个令人印象深刻的比较来自于开发mcp软件包：https://lindeloev.github.io/mcp/articles/packages.html的v博士。
下面，我使用您的样本时间序列，使用我自己开发的Rbeast软件包(显然是为了自我提升以及感知的亲切感)来生成一些快速的结果。Rbeast算法是一种Baysian变换点检测算法，它可以估计转换点发生的概率。它也可以用于将时间序列分解为季节性和趋势，但是很明显，您的时间序列是仅趋势的，因此在下面的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|
--------------------------------'Stack Overflow用户
发布于 2021-05-29 11:15:41
如果信号不太嘈杂，您可以使用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)
另一个似乎对您提供的数据更有效的选择是使用分段线性分割。
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
https://stackoverflow.com/questions/67749982
复制相似问题
腾讯云开发者
