Сглаживающий сплайн

Сглаживающий сплайн (англ. smoothing spline) это метод сглаживания (аппроксимации кривой набора зашумлённых исходных данных) с использованием сплайн-функций.

Определение

Пусть ( x i , Y i ) ; x 1 < x 2 < ⋯ < x n , i ∈ Z {displaystyle (x_{i},Y_{i});x_{1}<x_{2}<dots <x_{n},iin mathbb {Z} } последовательность наблюдений, порождённых выражением Y i = μ ( x i ) {displaystyle Y_{i}=mu (x_{i})} . Приближение сглаживающими сплайнами μ ^ {displaystyle {hat {mu }}} функции μ {displaystyle mu } определяется как функция (в классе дважды дифференцируемых функций), минимизирующая

∑ i = 1 n ( Y i − μ ^ ( x i ) ) 2 + λ ∫ x 1 x n μ ^ ″ ( x ) 2 d x . {displaystyle sum _{i=1}^{n}(Y_{i}-{hat {mu }}(x_{i}))^{2}+lambda int _{x_{1}}^{x_{n}}{hat {mu }}'(x)^{2},dx.}

Замечания:

λ ≥ 0 {displaystyle lambda geq 0} параметр сглаживания, контролирующий соотношение между точностью воспроизведения данных и «неровностью» аппроксимирующей функции.

интеграл вычисляется по всему диапазону x i {displaystyle x_{i}} .

при λ → 0 {displaystyle lambda o 0} (нет сглаживания), сглаживающий сплайн превращается в интерполяционный сплайн.

при λ → ∞ {displaystyle lambda o infty } (бесконечное сглаживание), штраф за неровность становится преобладающим и аппроксимация превращается в линейную МНК аппроксимацию.

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

в ранней литературе, с равноудалёнными x i {displaystyle x_{i}} , для вычисления штрафа вместо производной использовались конечные разности второго и третьего порядка.

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

Вывод сглаживающего сплайна

Для удобства разделим нахождение выражений, описывающих сглаживающий сплайн, на два этапа

Для начала найдём значения μ ^ ( x i ) ; i = 1 , … , n {displaystyle {hat {mu }}(x_{i});i=1,ldots ,n} .

Из этих значений найдём μ ^ ( x ) {displaystyle {hat {mu }}(x)} для всех x.

Начнём со второго шага.

Дан вектор m ^ = ( μ ^ ( x 1 ) , … , μ ^ ( x n ) ) T {displaystyle {hat {m}}=({hat {mu }}(x_{1}),ldots ,{hat {mu }}(x_{n}))^{T}} «подогнанных» значений, сумма квадратов в критерии сплайна — константа. Требуется только минимизировать ∫ μ ^ ″ ( x ) 2 d x {displaystyle int {hat {mu }}'(x)^{2},dx} , и минимизация — натуральный кубический сплайн, интерполирующий точки ( x i , μ ^ ( x i ) ) {displaystyle (x_{i},{hat {mu }}(x_{i}))} . Данный интерполяционный сплайн это линейный оператор, и он может быть представлен в форме:

μ ^ ( x ) = ∑ i = 1 n μ ^ ( x i ) f i ( x ) {displaystyle {hat {mu }}(x)=sum _{i=1}^{n}{hat {mu }}(x_{i})f_{i}(x)}

где f i ( x ) {displaystyle f_{i}(x)} набор базисных сплайн-функций. В результате штраф за негладкость имеет форму

∫ μ ^ ″ ( x ) 2 d x = m ^ T A m ^ . {displaystyle int {hat {mu }}'(x)^{2}dx={hat {m}}^{T}A{hat {m}}.}

где элементы A — ∫ f i ″ ( x ) f j ″ ( x ) d x {displaystyle int f_{i}'(x)f_{j}'(x)dx} . Базисные функции и матрица A зависят от конфигурации независимых переменных x i {displaystyle x_{i}} , но не от Y i {displaystyle Y_{i}} или m ^ {displaystyle {hat {m}}} .

Теперь, возвращаясь к первому шагу, взвешенная сумма квадратов может быть записана как

‖ Y − m ^ ‖ 2 + λ m ^ T A m ^ , {displaystyle |Y-{hat {m}}|^{2}+lambda {hat {m}}^{T}A{hat {m}},}

где Y = ( Y 1 , … , Y n ) T {displaystyle Y=(Y_{1},ldots ,Y_{n})^{T}} . минимизация по m ^ {displaystyle {hat {m}}} даёт

m ^ = ( I + λ A ) − 1 Y . {displaystyle {hat {m}}=(I+lambda A)^{-1}Y.}

Создание многомерных сплайнов

Из приведённого ограничения на формулу из определения x 1 < x 2 < ⋯ < x n {displaystyle x_{1}<x_{2}<dots <x_{n}} мы можем заключить, что алгоритм не работает для произвольного набора данных. Если планируется использование алгоритма для произвольного набора точек в многомерном пространстве необходим алгоритм, в котором нет таких ограничений. Возможное решение заключается во введении параметра таким образом, что входные данные могут быть представлены как одномерные функции, зависящие от данного параметра; затем можно применить сглаживание для каждой функции. В двумерном пространстве решение состоит в параметризации x {displaystyle x} и y {displaystyle y} как x ( t ) {displaystyle x(t)} and y ( t ) {displaystyle y(t)} где t 1 < t 2 < ⋯ < t n {displaystyle t_{1}<t_{2}<dots <t_{n}} . Подходящее решение для t {displaystyle t} это накопленное расстояние t i + 1 = t i + ( x i + 1 − x i ) 2 + ( y i + 1 − y i ) 2 {displaystyle t_{i+1}=t_{i}+{sqrt {(x_{i+1}-x_{i})^{2}+(y_{i+1}-y_{i})^{2}}}} где t 1 = 0 {displaystyle t_{1}=0} .

Более детальный анализ параметризации выполнен E.T.Y Lee.

Связанные методы

Сглаживающие сплайны имеют отношение, но отличаются от:

Регрессионных сплайнов. В этом методе данные апроксимируются набором базисных сплайн-функций с уменьшенным количеством узлов, как правило по МНК. Штрафы за негладкость не используются.
Penalized Splines. Сочетают уменьшенное количество узлов регрессионных сплайнов с штрафом за негладкость сглаживающих сплайнов.
Метод упругой карты. Данный метод сочетает МНК-штрафы для ошибки аппроксимации со штрафами за кривизну и растяжение аппроксимирующего множества и использует крупный шаг дискретизации для оптимизации проблемы.

Исходный код

Исходный код для сглаживающих сплайнов может быть взят из примеров к книге Carl de Boor’s A Practical Guide to Splines. Примеры написаны на Фортране. Обновлённые исходные коды также доступны на официальном сайте Carl de Boor’s [1].