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

О пакете BinSeqBstrap

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

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

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

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

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

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

2f22fc12ecde8e542a4fabcd261a11d5.png

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

n – общее количество точек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 отвергаетсяИтого — нулевая гипотеза на уровне значимости 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")

На удивление, все нашлосьНа удивление, все нашлось

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

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

© Habrahabr.ru