Регрессия и функции с неустранимыми разрывами первого рода

О пакете BinSeqBstrap

Постановка задачи

Допустим, у нас есть какая кусочно-гладкая функция f(x), к который прибавлен некий случайный шум, соответствующий условиям Гаусса-Маркова. И все хорошо, только эта самая функция f(x) – функция с неустранимым(-и)  разрывом(-ами) первого рода, то есть в какой-то точке левый и правый предел этой функции равны разным числам, а у функции есть скачок. Задача – как-то нужно научить алгоритм распознавать этот скачок.

Минутка теории

Теоретические основы изложены в виньетке, написанной Кэти МакДэйд и Флориана Пэйна из Кэмбриджа, опубликованной вот тут.

Если вкратце, то функционал пакета позволяет решить 4 разные подзадачи, базовая из которых - определение места и размера единственного скачка функции при известной ширине окна h

Основа: вот эта вот длинная формула

Результат этой формулы используется для вычисления места скачка

Сам метод основан на идеях непараметрической регрессии. Для определения местоположения точек разрыва выбирается окно из некоторого количества точек, и считается разница между суммами значений слева и справа. Для определения размера скачка используется разница значений сглаженной функции в определенной точек разрыва.

Практика

В коде это выглядит так:

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 - подогнанные значения

Часть 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")

Часть 3. Проверка статистической гипотезы о наличии скачка

Объявляем нулевую гипотезу о том, что размер скачка К равен 0 (альтернативная – не равен).

Дальше с помощью бутстрэпа получаем множество оценок величины этого скачка (много-много раз решаем чуть разные задачи из п.2) и смотрим, выполняется гипотеза или нет.

В коде это выглядит так:

test <- BstrapTest(y = y)
test$outcome # Результат проверки гипотезы о нулевом скачке
test$pValue

Часть 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")

Вместо заключения

Весьма интересная задумка. Для решения практических задач в текущем состоянии пакет применим мало, к тому же редко приходится данные аппроксимировать функциями с разрывами, но как инструмент вспомогательного анализа может пригодиться. Ну и интересно потестить, до каких пределов алгоритм работает хорошо, а когда начнет пропускать разрывы / находить лишние.