О пакете BinSeqBstrap
Постановка задачи
Допустим, у нас есть какая кусочно-гладкая функция f(x), к который прибавлен некий случайный шум, соответствующий условиям Гаусса-Маркова. И все хорошо, только эта самая функция f(x) – функция с неустранимым(-и) разрывом(-ами) первого рода, то есть в какой-то точке левый и правый предел этой функции равны разным числам, а у функции есть скачок. Задача – как-то нужно научить алгоритм распознавать этот скачок.
Минутка теории
Теоретические основы изложены в виньетке, написанной Кэти МакДэйд и Флориана Пэйна из Кэмбриджа, опубликованной вот тут.
Если вкратце, то функционал пакета позволяет решить 4 разные подзадачи, базовая из которых - определение места и размера единственного скачка функции при известной ширине окна h
Основа: вот эта вот длинная формула
![](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/2f2/2fc/12e/2f22fc12ecde8e542a4fabcd261a11d5.png)
Результат этой формулы используется для вычисления места скачка
![n – общее количество точек n – общее количество точек](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/659/71f/713/65971f7133171233706982528a243f66.png)
Сам метод основан на идеях непараметрической регрессии. Для определения местоположения точек разрыва выбирается окно из некоторого количества точек, и считается разница между суммами значений слева и справа. Для определения размера скачка используется разница значений сглаженной функции в определенной точек разрыва.
Практика
В коде это выглядит так:
library(BinSegBstrap)
## Часть 1
set.seed(1)
n <- 1:100
signal <- 0.1*n-5
signal[51:100] <- signal[51:100] + 5
y <- rnorm(n) + signal # Придумываем искусственные данные
est <- estimateSingleCp(y = y, bandwidth = 0.1)
est$cp # определяем точку разрыва
est$size # определяем скачок функции
plot(y, pch = 16, col = "grey30")
lines(signal)
lines(est$est, col = "red") # est$est - подогнанные значения
![Местоположение скачка точно установлено, его размер – примерно. Местоположение скачка точно установлено, его размер – примерно.](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/0b6/47a/b32/0b647ab32ed0bdf51936498811f9ca99.png)
![А вот графика. Красная линия - это сглаженная регрессия, черная - подобранные линии уравнений А вот графика. Красная линия - это сглаженная регрессия, черная - подобранные линии уравнений](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/ca8/723/7d6/ca87237d60fc67f1d6b218b3e3fe0a73.png)
Часть 2. Определение параметров скачка при неизвестной величине
Собственно, к исходной формуле достраивается только алгоритм кросс-валидации, который позволяет найти лучший размер окна, в котором искать скачок.
Зададим алгоритму задачку посложнее:
set.seed(1)
n <- 1:100
signal1 <- 5-0.1*n
signal2 <- 0.01*n^2-20
signal<-c(signal1[1:50],signal2[51:100])
y <- rnorm(n) + signal
est <- estimateSingleCp(y = y)
est$bandwidth # подобранная ширина окна
est$cp # определяем точку разрыва
est$size # определяем скачок функции
plot(y, pch = 16, col = "grey30") # plot of observations, true and estimated signal
lines(signal)
lines(est$est, col = "red")
![Снова все определено верно Снова все определено верно](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/5df/4b3/c78/5df4b3c7890b3a5ce3775ea19a9cfb6f.png)
Часть 3. Проверка статистической гипотезы о наличии скачка
Объявляем нулевую гипотезу о том, что размер скачка К равен 0 (альтернативная – не равен).
Дальше с помощью бутстрэпа получаем множество оценок величины этого скачка (много-много раз решаем чуть разные задачи из п.2) и смотрим, выполняется гипотеза или нет.
В коде это выглядит так:
test <- BstrapTest(y = y)
test$outcome # Результат проверки гипотезы о нулевом скачке
test$pValue
![Итого – нулевая гипотеза на уровне значимости 0.05 отвергается Итого – нулевая гипотеза на уровне значимости 0.05 отвергается](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/a20/f8c/358/a20f8c3585b7cacc00b2e31ba30d2388.png)
Часть 4. Поиск всех возможных точек разрыва
Составим функцию, в которой есть две точки разрыва, и посмотрим, справится ли алгоритм с ней
set.seed(1)
signal1<-abs(-8-seq(from=-10, to=-5, length.out=50))
signal2<-3*exp(seq(from=-5, to=0, length.out=50))
signal3<-sqrt(seq(from=0, to=10, length.out=50))-4
signal <- c(signal1,signal2,signal3)
y <- rnorm(150) + signal
est <- BinSegBstrap(y = y)
est$cps # находим точки разрыва
plot(y, pch = 16, col = "grey30") # plot of observations, true and estimated signal
lines(signal)
lines(est$est, col = "red")
![На удивление, все нашлось На удивление, все нашлось](https://webcf.waybackmachine.org/web/20220305121122/https://habrastorage.org/getpro/habr/upload_files/790/94d/713/79094d713394d2b95ceabdc7cf7aba88.png)
Вместо заключения
Весьма интересная задумка. Для решения практических задач в текущем состоянии пакет применим мало, к тому же редко приходится данные аппроксимировать функциями с разрывами, но как инструмент вспомогательного анализа может пригодиться. Ну и интересно потестить, до каких пределов алгоритм работает хорошо, а когда начнет пропускать разрывы / находить лишние.