обнаружение изменения точки накопительного распределения

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-linear-regression-model

но в R. нет эквивалентного пакета.

Кто-нибудь может подсказать, как это сделать?

Вот ожидаемый результат.

Answer 1

Это не ответ, а комментарий (слишком длинный, чтобы его можно было редактировать в разделах комментариев).

Я считаю ваш числовой пример интересным, особенно для сравнения с результатами, полученными благодаря методу изстатья: 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 2

Мы можем использовать пакет 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")
   

Более подробную информацию можно найти в справочном руководстве segmented на CRAN.Возвращаемый объект fit.seg также содержит оценки параметров для каждого куска.