CC BY
66
УДК 519.234
Об одном методе сглаживания двумерной поверхности
П. Г. Любин
ФБГОУ ВО МГТУ «СТАНКИН», Москва,, Россия
Регрессионный анализ ставит перед собой задачу отыскания функциональной зависимости между наблюдаемыми величинами изучаемого процесса. При этом исходные данные являются реализацией случайной величины, поэтому рассматривается зависимость математического ожидания. Такую задачу можно решать путём «сглаживания» исходных данных. Под сглаживанием понимается попытка удаления шума и несущественных фрагментов при сохранении наиболее важных свойств структуры данных, то есть результат подобен математическому ожиданию. Сглаживание данных, как правило, осуществляется путём параметрической или непараметрической регрессии. В случае параметрической регрессии необходимы априорные знания о форме уравнения регрессии. Большинство исследуемых данных, однако, невозможно параметризовать. С этой точки зрения непараметрическая и полупараметрическая регрессии представляются лучшим подходом к решению задачи сглаживания. Целью исследования ставилось разработка и реализация алгоритма быстрого сглаживания двумерных данных. Для достижения этой цели были проанализированы предыдущие работы в данной области и разработан свой подход, улучшающий предыдущие. В результате, в данной работе представлен алгоритм, который быстро и с минимальным потреблением памяти очищает данные от «шума» и «несущественных» частей. Для подтверждения «эффективности» алгоритма проведены сравнения с другими общепризнанными подходами на смоделированных и реальных данных. Результаты этих сравнений также приведены в статье.
Ключевые слова: непараметрическая регрессия, двумерное сглаживание, штрафные сплайны, сглаживающие сплайны, скользящий контроль, двумерное дискретное косинусное преобразование
При анализе некоторого реального процесса исходные данные зашумлены, из-за чего необходимо «сглаживание». Под сглаживанием понимается попытка фильтрации шума или несущественных фрагментов при сохранении наиболее важных свойств структуры данных. Рассмотрим следующую модель
где £ — гауссов белый шум. Предполагается, что функция f (х) должна быть гладкой, т.е. иметь непрерывные производные до некоторого порядка. Сглаживание данных, как правило, осуществляется путём параметрической или непараметрической регрессии. В случае параметрической регрессии необходимы некоторые априорные знания о форме уравнения регрессии, которое достаточно хорошо описывало бы исходные данные. Большинство наблюдаемых данных, однако, невозможно параметризовать с точки зрения задания функции / (х) в аналитическом виде. С этой точки зрения непараметрическая и полупараметрическая регрессии являются лучшим подходом к решению задачи (1). Одним из классических подходов к сглаживанию является использование различных модификаций метода наименьших квадратов со штрафом. Он впервые был продемонстрирован вначале 1920-х в работе [1] и подробно разобран в книге [2]. Эта техника заключается в минимизации некоторого функционала, который уравновешивает «адекватность» и «гладкость» оценки
1. Постановка задачи
У = У + £,
(1)
^(у) = Явв + X ■ Р(у) = ||у - у||2 + Л ■ Р(у),
(2)
где || ■ || — евклидова норма. Параметр Л является вещественным положительным числом, контролирующим гладкость решения: при его возрастании гладкость у также растёт. Когда штрафная функция записана в виде интеграла квадрата производной р-порядка, регрессия называется сглаживающим сплайном [1,3,4]. Другим простым и эффективным подходом к решению задачи (1 )является квадратичный вид штрафной функции в следующей форме [5]:
р (у) = \\т2, (3)
где И — трёхдиагональная матрица следующего вида
--1 1 1 -2 1
1 -2 1 1 -1
Данная статья является продолжением исследований автора, затронутых в статьях [6,7], а также некоторые идеи почерпнуты из статьи [8].
2. Одномерное сглаживание
Пусть имеются равноотстоящие точки {жг}к<п и значения функции отклика в этих точках следующего вида
Уi = / (Хг) + £г, (4)
где £г ~ N(0,а2). Пусть у является оценкой функции /(Хг). Тогда минимизация выражения (2) приводит к следующему уравнению
у = Н(Л) ■ у,
(5)
где Н(X) = (I + Л ■ БТБ)-1 представляет собой проекционную матрицу, а Л — это сглаживающий параметр. Сглаживающий параметр выбирается путём минимизации функционала следующего вида
ССУ (X) =
ВВв (Х)/п
(1 - Тг(Я(Х))/п)2'
(6)
Такой подход называется методом скользящего контроля (сго88уаНёа1;юп). При равноотстоящих наблюдениях можно использовать свойства матрицы И, с помощью которых можно упростить вычисления ОСУ. Это свойство заключается в возможности разложить матрицу И в следующее произведение иГит, где матрица и является унитарной и представляет собой дискретное косинусное преобразование [9]. Тогда Явв можно переписать в следующем виде:
КББ = ||у - у\\2 = \\Н(X) ■ у - у\\2 = || ((/ + Л ■ БТО)- - /) ■ у
1
И' + Л-Г2)- - 7) ■ ^ = £ - 1)
- 0 ОСТ^(р).
2
2
В таком случае функционал (6) принимает более удобный для вычислений вид:
п ■ -да - 1) -пстг2(у) ОСУ(X) =--, \ ^-. (7)
п
3. Двумерное сглаживание
Пусть имеются равномерная сетка {(х\¿,х2^)}т«П1 1<з<П2 и значения функции отклика в узлах этой сетки следующего вида
У1,3 =1(Х\) + £1,3, (8)
где £i,j ~ N(0,а2). В таком случае можно представить значения функции отклика в виде матрицы У, в которой у^^ является элементом г-й строки j-го столбца. Тогда
сглаженные значения будем обозначать как У. Введём операцию нес, которая заключается в записи матрицы в виде вектор-столбца путём выкладывания столбцов матрицы друг за другом. Очевидно, что оценка Vес(У) имеет следующий вид:
Vес(У) = ( Н Х2 ®НХ1) ■ Vес(У) = НХ2уХ1 ■ Vес(У), (9)
где НХ1, НХ2 — проекционные матрицы соответствующего измерения. Очевидно, что эти проекционные матрицы имеют следующий вид
НХг = (1Пг + ХгР>тпрп% )-1 , г = 1, 2. (10)
Применяя вышеизложенный подход и свойства тензорного произведения [10], выражение (9) можно упростить следующим образом
у=(НХ2 ®НХ1) -у =(1П2 + Х20120П2)-1 ® (1П1 +Х1В11БП1 )-1 -у =
= "•' {. ттЬг. И (ттк И ■» =
= иХ2 ® иХ1 ■ ® ( 1 + \112 ) ■ и*2 ® и1 ■ у = ■ ■ и12,х1 ■ у.
Для автоматического поиска оптимальных значений А1 и Л2 предлагаем также
использовать ОСУ адаптированный для двумерного случая:
ОСУ = ■ (11)
Принимая во внимание свойство следа тензорного произведения матриц [10], получаем Тг(НХ2 Х1) = ^2 т + X 2 ■ Е -| IX1 2 . Очевидно, что основным затратным местом в части вычислений является расчёт Явв, так как требуется вычисление у для всех комбинаций А1 и А2. Данный расчёт можно упростить следующим образом:
КББ = ||у - у\\2 = \\НХ2,Х1 ■ у - у\\2 = \\(НХ2,Х1 - 1п) ■ у\\2 =
■ (Г Х2,Х1 - 1п) ■ их2,х1 ■ у\\ = (и^,х1 ■ (ГХ2,Х1 - 1п) ■ их2,х1 ■ у) ■ {их2,Х1 ■ (Гх2,Х1 - 1п) ■ их2,х1 ■ у) =
= (БСТ2 ■ у)Т ■ (ГЖ2)Л1 - 1п)2 ■ БСТ2 ■у = ^ (7х2,Ж1 - 1)2 ■ (ВСТ2 ■ у)2 .
Здесь ИСТ2 — это двумерный аналог дискретного косинусного преобразования. Из упрощённой формулы видно, что преобразование необходимо выполнить один раз, а меняются только значения 7Х2,Х1 в зависимости значений А1 и А2. Данный подход реализован на языке И. Для демонстрации преимуществ изложенного подхода выполнены численные эксперименты: на искусственно сгенерированных данных и на реальных данных.
4. Численный эксперимент
Для проведения эксперимента смоделирована следующая задача: взята функция 8ш(2-^(ж-0, 5)3) ■ ео8(4^у) и зашумлена случайными значениями из нормального распределения N(0, 0, 22) (рис. 1). Сглаживание проводилось изложенным подходом и при помощи пакета МССУ [11], в котором реализовано сглаживание штрафными сплайнами, в том числе и для многомерного случая с использованием тензорного произведения базовых функций. Ниже приведена таблица с результатами сглаживания на различных сетках.
Рис. 1. Функция 8ш(2-^(ж — 0, 5)3) • 008(4"^) с шумом
Таблица 1
Результаты сглаживания смоделированной задачи на разных сетках
Параметр Р-сплайны с ДКП САМ с 102 узлами САМ с 202 узлами
РЯБ МЯЕ Корреляция с истин. значениями Время оцен. (с) 9, 488243 0, 001483 0, 9993394 1, 941 11,72485 0, 001832 0, 996919 10,237 9, 87163 0, 00154 0, 9991624 29, 875
Для демонстрации практического применения вышеизложенного подхода взяты данные о смертности в России. Данные взяты из открытого источника [12] и содержат наблюдения за 1959-2010 годы по возрастам 0-110. Оценивание проводилось на части данных, которая относится к старшим возрастам (50-101, рис. 2). К этим годам выборка достаточно сильно уменьшается, в результате чего наблюдения содержат ошибки и выбросы, из-за которых не видно общей картины происходящих процессов. Таким образом, анализируемые данные представляют собой равномерно расположенные значения коэффициентов смертности на сетке размерности 52 х 52. Сглаживание проводилось изложенным подходом, пакетом МССУ и параметрической моделью Ли-Картера, которая стала классической при оценивании двумерной поверхности смертности. Далее приведена таблица с результатами оценивания.
Рис. 2. Коэффициенты смертности населения России старше 50 лет с 1959 г.
по 2010 г.
Таблица 2
Результаты сглаживания двумерной поверхности смертности России
Параметр Р-сплайны с ДКП Модель Ли-Картера САМ с 122 узлами
МЯЕ Время оцен. (с) 0,21637 0, 0000905 0, 49 18,5092 0, 0077379 1,194 0, 41395 0, 0001731 4,185
5. Заключение
Из приведённых результатов очевидно, что подход описанный в данной статье более эффективен, так как имеет большую скорость сглаживания и малое потребление памяти. Также хочется отметить, что при росте выборки скорость оценивания растёт незначительно при сохранении качества оценки. Результатами проделанной работы являются:
1. получены выражения для двумерного случая с учётом двух параметров сглаживания;
2. полученный подход реализован на языке И;
3. выполнено сравнение подхода с другими аналогичными подходами и моделями.
В последующих работах предполагается рассмотрение следующих возможностей:
— расширение метода на многомерный случай;
— использование других распространённых критериев выбора сглаживающих параметров, например, В1С или А1С;
— использование более быстрого метода поиска минимального значения функционала ОСУ.
Литература
1. Whittaker E. T. On a New Method of Graduation // Proceedings of the Edinburgh Mathematical Society. — 1923. — Vol. 41. — Pp. 62-75.
2. Wahba G. Spline Models for Observational Data. — Society for Industrial Mathematics, 1990. — ISBN 9780898712445.
3. Schoenberg I. J. Spline Functions and the Problem of Graduation // Proceedings of the National Academy of Sciences of the United States of America. — 1964. — Vol. 52. — Pp. 947-950.
4. Takezawa K. Introduction to Nonparametric Regression. — Wiley & Sons, Inc., 2005. — ISBN 9780471745839.
5. Weinert H. L. Efficient Computation for Whittaker-Henderson Smoothing // Computational Statistics & Data Analysis. — 2007. — Vol. 52. — Pp. 959-974.
6. Щетинин Е. Ю., Любин П. Г. Робастный алгоритм построения сглаживающих сплайнов // Научное обозрение. — 2015. — Т. 1. — С. 86-94.
7. Любин П. Г., Щетинин Е. Ю. Стохастические модели сглаживания и прогнозирования коэффициентов смертности // Научное обозрение. — 2015. — Т. 18. — С. 147-155.
8. Xiao L., Li Y., Ruppert D. Fast Bivariate P-splines: the Sandwich Smoother // Journal of the Royal Statistical Society. — 2013. — Vol. 75. — P. 577-599.
9. Garcia D. Robust Smoothing of Gridded Data in One and Higher Dimensions with Missing Values // Computational Statistics & Data Analysis. — 2010. — Vol. 54. — P. 1167-1178.
10. Seber G. A Matrix Handbook for Statisticians. — Wiley-Interscience, 2007. — ISBN 9780471748694.
11. Wood S. MGCV: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness Estimation. — R package version 1.8.10. https://cran.r-project.org/web/ packages/mgcv/index.html, r package version 1.8.10.
12. The Human Mortality Database. — http://www.mortality.org/.
UDC 519.234
On a Method of Two-Dimensional Smoothing P. G. Lyubin
Moscow State Technology University "STANKIN", Moscow, Russia
Regression analysis has the task of finding a functional relationship between the observed values the studied process. The raw data is the realization of a random variable, it is therefore considered dependent on the expectation. This problem can be solved by "smoothing" the raw data. Smoothing is the process of removing the noise and insignificant fragments while preserving the most important properties of the data structure. It is similar to finding the expectation of data. Data smoothing usually attained by parametric and nonparametric regression. The nonparametric regression requires a prior knowledge of the regression equation form. However,
most of the investigated data cannot be parameterized simply. From this point of view, non-parametric and semiparametric regression represents the best approach to smoothing data. The aim of the research is development and implementation of the fast smoothing algorithm of two-dimensional data. To achieve this aim previous works in this area have been analyzed and its own approach has been developed, improving the previous ones. As a result, this paper presents the algorithm that quickly and with minimal memory consumption cleanses the data from the "noise" and "insignificant" parts. To confirm the "efficiency" of the algorithm the comparisons with other generally accepted approaches were carried out on simulated and real data with other generally accepted approaches. The results of these comparisons are also shown in the paper.
Key words and phrases: nonparametric regression, two-dimensional estimation, penalized splines, smoothing splines, cross-validation, discrete cosine transform
References
1. E. T. Whittaker, On a New Method of Graduation, Proceedings of the Edinburgh Mathematical Society 41 (1923) 62-75.
2. G. Wahba, Spline Models for Observational Data, Society for Industrial Mathematics, 1990.
3. I. J. Schoenberg, Spline Functions and the Problem of Graduation, Proceedings of the National Academy of Sciences of the United States of America 52 (1964) 947-950.
4. K. Takezawa, Introduction to Nonparametric Regression, Wiley & Sons, Inc., 2005.
5. H. L. Weinert, Efficient Computation for Whittaker-Henderson Smoothing, Computational Statistics & Data Analysis 52 (2007) 959-974.
6. E. Y. Shetinin, P. G. Lyubin, Robust Smoothing with Splines, Science Review 1 (2015) 86-94, in Russian.
7. P. G. Lyubin, E. Y. Shetinin, Stochastic Models of Mortality Estimation, Science Review 18 (2015) 147-155, in Russian.
8. L. Xiao, Y. Li, D. Ruppert, Fast Bivariate P-splines: the Sandwich Smoother, Journal of the Royal Statistical Society 75 (2013) 577-599.
9. D. Garcia, Robust Smoothing of Gridded Data in One and Higher Dimensions with Missing Values, Computational Statistics & Data Analysis 54 (2010) 1167-1178.
10. G. Seber, A Matrix Handbook for Statisticians, Wiley-Interscience, 2007.
11. S. Wood, MGCV: Mixed GAM Computation Vehicle with GCV/AIC/REML Smoothness Estimation, r package version 1.8.10.
URL https://cran.r-project.org/web/packages/mgcv/index.html
12. The human mortality database, last visited on 25.02.2016. URL http://www.mortality.org/
© Любин П.Г., 2016
