Регрессия и функции с неустранимыми разрывами первого рода
О пакете BinSeqBstrap
Постановка задачи
Допустим, у нас есть какая кусочно-гладкая функция f (x), к который прибавлен некий случайный шум, соответствующий условиям Гаусса-Маркова. И все хорошо, только эта самая функция f (x) — функция с неустранимым (-и) разрывом (-ами) первого рода, то есть в какой-то точке левый и правый предел этой функции равны разным числам, а у функции есть скачок. Задача — как-то нужно научить алгоритм распознавать этот скачок.
Минутка теории
Теоретические основы изложены в виньетке, написанной Кэти МакДэйд и Флориана Пэйна из Кэмбриджа, опубликованной вот тут.
Если вкратце, то функционал пакета позволяет решить 4 разные подзадачи, базовая из которых — определение места и размера единственного скачка функции при известной ширине окна h
Основа: вот эта вот длинная формула
Результат этой формулы используется для вычисления места скачка
n — общее количество точек
Сам метод основан на идеях непараметрической регрессии. Для определения местоположения точек разрыва выбирается окно из некоторого количества точек, и считается разница между суммами значений слева и справа. Для определения размера скачка используется разница значений сглаженной функции в определенной точек разрыва.
Практика
В коде это выглядит так:
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
Итого — нулевая гипотеза на уровне значимости 0.05 отвергается
Часть 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")
На удивление, все нашлось
Вместо заключения
Весьма интересная задумка. Для решения практических задач в текущем состоянии пакет применим мало, к тому же редко приходится данные аппроксимировать функциями с разрывами, но как инструмент вспомогательного анализа может пригодиться. Ну и интересно потестить, до каких пределов алгоритм работает хорошо, а когда начнет пропускать разрывы / находить лишние.
