累积分布的变化点检测

Question

我有一个累积降雨量时间序列，我想检测变化点。这是数据。

structure(list(DAY = 1:365, CUMSUM = c(0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.8, 6.9, 6.9, 6.9, 
6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 7.4, 
7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 7.4, 22.6, 
22.6, 22.6, 22.6, 22.6, 22.6, 22.8, 26.7, 41.3, 41.3, 44.7, 44.7, 
44.7, 86.8, 92.6, 92.6, 115.2, 117, 126, 134.9, 134.9, 134.9, 
140.7, 140.7, 140.7, 146.5, 146.7, 146.7, 151.7, 152.7, 196.5, 
242.7, 293.4, 331.4, 340, 345.6, 369.5, 442.6, 459, 464.6, 464.6, 
468.2, 475.6, 484.2, 487.8, 487.8, 511, 515, 515, 515, 528.8, 
547.6, 549.4, 549.8, 550, 552.4, 585.9, 798.5, 1062.5, 1107.9, 
1124.5, 1154, 1169.4, 1416.4, 1457.6, 1457.6, 1457.6, 1461.2, 
1464, 1524.7, 1539.5, 1552, 1592.8, 1599.4, 1608.6, 1611.6, 1616.2, 
1656.6, 1667.6, 1667.6, 1668.8, 1680, 1687.1, 1697.9, 1704.7, 
1726.6, 1726.6, 1727.6, 1732.6, 1750.2, 1834.4, 1882.2, 1915.6, 
1940, 1976.6, 2001.2, 2026.4, 2042.6, 2078.1, 2101.2, 2109.2, 
2109.2, 2109.2, 2109.2, 2117, 2117, 2120.2, 2142.4, 2153.4, 2173.4, 
2174.4, 2174.4, 2174.4, 2178.4, 2213.5, 2365.1, 2449.7, 2565.5, 
2673.7, 2749.9, 2830.3, 2896.2, 2920.8, 3236.4, 3266.8, 3288.9, 
3371.5, 3428.5, 3642.5, 3764.9, 3774.9, 3818.7, 3818.7, 3830.9, 
3953.7, 4127.8, 4206, 4217.7, 4217.7, 4219.9, 4220.9, 4220.9, 
4361.1, 4378, 4378, 4388.4, 4393.4, 4417.3, 4419.9, 4419.9, 4419.9, 
4470.3, 4480.3, 4480.7, 4490.7, 4492.9, 4493.4, 4504, 4504, 4504, 
4505.4, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 4509.8, 
4510.4, 4510.4, 4512.8, 4515.4, 4517.8, 4527.5, 4532.1, 4539.7, 
4541.7, 4573.3, 4606.5, 4607.3, 4613.5, 4613.5, 4613.5, 4613.5, 
4613.5, 4613.5, 4613.5, 4613.5, 4613.5, 4613.5, 4613.9, 4621.1, 
4621.1, 4621.1, 4636.5, 4647.9, 4649.1, 4649.3, 4649.3, 4649.3, 
4655, 4655, 4663.6, 4663.6, 4664.2, 4664.2, 4665, 4665, 4665, 
4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 4665, 
4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 
4665.9, 4665.9, 4665.9, 4665.9, 4665.9, 4673.1, 4673.1, 4673.1, 
4673.1, 4673.1, 4673.1, 4673.1, 4673.1, 4673.1, 4673.5, 4673.5, 
4673.5, 4673.5, 4673.5, 4673.5, 4673.5, 4673.5, 4673.5)), .Names = 
c("DAY","CUMSUM"), class = "data.frame", row.names = c(NA, -365L))

我想应用两相线性回归来使用R检测变化点。

这里有一个matlab代码 https://www.mathworks.com/matlabcentral/fileexchange/26804-two-phase-线性-回归-模型

但 R 中没有等效的包。

有人可以建议如何做到这一点吗？

这是预期的输出。

Answer 1

我们可以使用R包

segmented

；这是一个分步示例。

1. 加载库。

   library(segmented);
   

2. 将具有两个断点的分段线性模型与样本数据拟合（这里我假设

   df

   包含数据作为

   data.frame

   ）。请注意，我们必须对断点提供一些猜测。

   fit <- lm(CUMSUM ~ DAY, data = df);
   fit.seg <- segmented(fit, psi = c(100, 200));
   fit.seg;
   #Call: segmented.lm(obj = fit, psi = c(100, 200))
   #
   #Meaningful coefficients of the linear terms:
   #(Intercept)          DAY       U1.DAY       U2.DAY
   #     -58.20         1.25        35.70       -34.98
   #
   #Estimated Break-Point(s):
   #psi1.DAY  psi2.DAY
   #   153.8     272.9
   

3. 我们绘制曲线并用红色标记断点估计。

   library(ggplot2);
   ggplot(df, aes(DAY, CUMSUM)) +
       geom_line() +
       geom_vline(data = as.data.frame(fit.seg$psi), aes(xintercept = `Est.`), col = "red")
   

4. 更多详细信息可以在 CRAN 上的

   segmented

   参考手册中找到。返回对象

   fit.seg

   还包含每个部分的参数估计。

Answer 2

这不是答案，而是评论（太长，无法在评论部分进行编辑）。

我发现你的数值示例很有趣，特别是与通过论文方法获得的结果进行比较：https://fr.scribd.com/document/380941024/Regression-par-morceaux-Piecewise-Regression-pdf

第 30-31 页给出的算法不是迭代的，不需要初始猜测。结果如下第一张图所示：

拟合的分段函数由三个线性段组成。但第一段和第三段并不完全按照您的问题中的要求水平。

事实上，这来自于参考论文中提到的积分方程的拟合。为了获得水平第一和第三段，必须使用参数 p1=p3=0 来简化微积分。此外，参数q1=0和q3=4673.5是先验已知的。算法得到简化：

结果是：

结果与 R 包的结果略有不同：a1=153.8 和 a2=272.9

有趣的是，最接近的结果是假设第一段和第三段不完全水平（a1=152 和 a2=274）。

当然，获得略有不同的结果并不奇怪，因为在每种情况下，回归的标准都不相同（并且我们不确切知道它们在 R 包中是什么）。

Answer 3

我们可以从不同的角度来看待问题。

通过比较数据，我们得到了图：。

我们可以使用一个简单的均值变化模型：

> fastcpd::fastcpd.mean(diff(data))
                                                                                                                                      
A fastcpd object.
Available methods to evaluate the object are:
plot, print, show, summary


Change points:
[1] 173 242

或者，也可以在此处安装 AR(1) 模型：

> suppressMessages(fastcpd::fastcpd.ar(diff(data), 1))
                                                                                                                                      
A fastcpd object.
Available methods to evaluate the object are:
plot, print, show, summary


Change points:
[1] 173 236

免责声明：我是 R 包的开发者

fastcpd

。