代码之家 › 专栏 › 技术社区 › Lyndz

累积分布的变点检测

cdf linear-regression r

Lyndz · 技术社区 · 6 年前

我有一个累积的降雨时间序列,我想检测变化点。以下是数据。

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,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.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,
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,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),姓名=
c(“day”,“cumsum”),class=“data.frame”,row.names=c(na,-365l))
我想在这里使用r检测变化点时应用一个两相线性回归。


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

但是在R中没有同等的包装。

有人能建议怎么做吗?

这是预期的输出。
.


我想申请两相线性回归在这里用R检测变化点。

这里有一个matlab代码
https://www.mathworks.com/matlabcentral/fileexchange/26804-two-phase-linear-regression-model

但是在R中没有等价的包。

有人能建议怎么做吗?

这是预期的输出。

2 回复 | 直到 6 年前

Maurits Evers 6 年前

库(分段);

fit<-lm(cumsum~day,data=df);
#(截获)第u1天第u2天
#

几何线条()。+




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



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


   
    
     
    
   
   
   
   
    
     
      分段
     
     
      reference manual
     
     
      fit.seg

JJacquelin 6 年前

这不是答案,而是注释(太长,无法在“注释”部分中编辑)。

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

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

结果是:

https://fr.scribd.com/document/380941024/Regression-par-morceaux-Piecewise-Regression-pdf

事实上,这是从一个积分方程的拟合,如文中提到的。为了得到水平的第一段和第三段,必须用参数p1=p3=0简化微积分。此外,参数q1=0和q3=4673.5是先验的。算法简化:

结果与R包稍有不同:A1=153.8,A2=272.9。

当然,获得稍微不同的结果并不奇怪,因为在每种情况下,回归的标准都不相同(我们不知道它们在R包中是什么)。