ВЫБОР ПАРАМЕТРА СГЛАЖИВАНИЯ КУБИЧЕСКОГО СПЛАЙНА МЕТОДОМ ПЕРЕКРЕСТНОЙ ПРОВЕРКИ Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

MATHEMATICAL SCIENCES

CHOISE OF THE CUBIC SPLINE SMOOTHING PARAMETER BY THE CROSS-VALIDATION

METHOD

Voskoboinikov Yu.

Doctor of Physical and Mathematical Sciences, Professor, Head of the Department of Applied Mathematics Novosibirsk State University of Architecture and Civil Engineering (Sibstrin), Novosibirsk

ВЫБОР ПАРАМЕТРА СГЛАЖИВАНИЯ КУБИЧЕСКОГО СПЛАЙНА МЕТОДОМ

ПЕРЕКРЕСТНОЙ ПРОВЕРКИ

Воскобойников Ю.

доктор физико-математических наук, профессор, заведующий кафедрой прикладной математики Новосибирского государственного архитектурно-строительного

университета (Сибстрин)

Abstract

The smoothing cubic spline of a defect unit is a fairly versatile tool for filtering measurement noise. The filtering (smoothing) error depends significantly on the value of the so-called smoothing parameter. In this paper, an algorithm for estimating the optimal parameter based on the cross-validation method is constructed and investigated in detail. In this case, special attention is paid to the estimation of the optimal parameter in conditions of correlated noise.

Аннотация

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

Keywords: noise filtering, smoothing cubic spline, smoothing parameter, parameter choise by cross-validation method.

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

Введение и задачи исследования. Сглаживание и дифференцирование зашумленных сигналов (данных) являются наиболее распространёнными задачами при обработке экспериментальных данных [1], в частности при решении задач непараметрической идентификации [2-4]. Следует напомнить, что задача дифференцирования является некорректно поставленной задачей [5]. Для фильтрации шума измерений и устойчивого дифференцирования часто используют сглаживающие кубические сплайны (СКС) дефекта единица. СКС имеет на всем интервале определения сплайна непрерывные первую и вторую производные. Главной проблемой при использовании СКС на практике является выбор параметра сглаживания, от величины которого существенно зависит ошибки сглаживания. Очевидно желание выбрать такое значение (назовем его оптимальным), при котором величина ошибки была бы минимальной. Следует выделить два случая: когда дисперсия шума измерения достоверно известна, и когда такая информа-

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

Поэтому целью данной работы является построение и исследование алгоритма выбора параметра сглаживания при неизвестной дисперсии шума измерения и коррелированном шуме измерения на основе метода перекрестной проверки. При этом решаются две задачи:

• построение алгоритма выбора параметра сглаживания сглаживающего кубического сплайна при неизвестной дисперсии;

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

Сглаживающий кубический сплайн. Сначала приведем некоторые определения и понятия сглаживающего кубического сплайна (СКС), необходимые для дальнейшего изложения алгоритма выбора параметра сглаживания СКС [6,7].

Предположим, что на некотором интервале

[7, 7 ] заданы п узлов 7] = ^ < ¿2 <... < ¿п = 7 , и в этих узлах измерены значения функции (сигнала) / (¿) :

fi=f(ti) + Л/, =

(1)

где ^ - случайный шум измерений с нулевым

средним и дисперсией а^ (равноточные измерения). Функция Би (Ц) называется кубическим сплайном дефекта единица, если:

а) на каждом отрезке , +1) функция Бп (Ц) является кубическим многочленом вида:

Бп () = а; + Ь ( - ^ ) + С; ( - Ц )2+ ё; ( - Ц )3;(2)

б) функция Бп (Ц) дважды непрерывно дифференцируема на всем интервале [7 ,72 ].

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

кубический сплайн Бп а (Ц) в общем случае не проходит через эти точки, а проходит в некоторой окрестности (ее размеры определяются величиной параметра сглаживания а) от этих точек, обеспечивая тем самым сглаживание (фильтрацию) шума измерений. Для следующих часто используемых на практике краевых условий [6,7]:

• условия на значения первой производной

вида

Я'п,а (Ч ) = 51, ,а (п ) = ¿п ; (3)

• условия на нулевые значения второй производной (так называемые естественные краевые условия)

= 0, ЗД ) = 0, (4)

СКС доставляет минимум функционалу [6,7]:

*п п

(5) = а{ (0|2 Л + X Р~]1 (У - % ))2 .(5)

где р; - весовые множители, отражающие

точность /-го измерения (в случае равноточных измерений задаются одинаковыми, например,

Р1 = 1).

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

параметра фильтрация шумов незначительна, и вычисленная первая производная имеет осциллирующую форму, что характерно для неустойчивых решений некорректно поставленных задач [5]. При больших значениях параметра график сплайна стремится к прямой линии (вторая производная стремится к нулю), и происходит существенное «сглаживание» информативных составляющих обрабатываемой функциональной зависимости. Между этими двумя предельными значениями находится значение параметра (назовем его оптимальным), для которого ошибка сглаживания (в принятой норме) минимальна. К сожалению, точное значение оптимального параметра на практике вычислить не удается (из-за отсутствия априорной информации о значениях или характеристиках «точного» сигнала). Поэтому попытаемся построить оценку для оптимального параметра сглаживания, адаптируя для этой цели метод перекрестной проверки (или метод перекрестной значимости), в зарубежной литературе - cross-validation method (CVM). Метод использовался при построении моделей (например, [8]) и для выбора параметра регуляризации в алгоритмах решения некорректных задач (например, [9, стр. 77-79]). Идейная основа этого метода видна из следующих рассуждений [1, стр. 68-69].

Оценивание оптимального параметра сглаживания. Предположим, что имеются два вектора

У и У размерности п. составленные из двух

наборов наблюдений (2). По первому вектору f построен сплайн Sn а (t) , из значений которого в узлах сформирован вектор sn а . Введем функцио-

нал:

U (а)= M

2* f ~S:

га,а

(6)

где математическое ожидание М [•] берется по

ансамблю распределения шума измерений ^ , а ||| означает евклидову норму вектора. Можно показать, что и (а) допускает представление:

[/(а)=А/[||/*-/||2 + 2(/*-/)Г.(/-^а)+||/-^а||2]Х7)

где / - вектор, составленный из «точных» значений У = У(7,-), / = 1.../7. Если матсматичс-

ское ожидание шума измерений равно нулю и шумы измерений двух набора данных не коррели-рованы между собой, то (7) можно переписать как:

U (а) = no2+ M где второе слагаемое (а

- s

n, а

= па,

,+Д(а) ,(8)

Д(а) = M

- s„

(9)

определим, как среднеквадратическую ошибку (СКО) сглаживания. Так как первое слагаемое в (8)

не зависит от параметра сглаживания а, то минимум и (а) соответствует минимуму СКО, т.е. задача вычисления оптимального параметра ао^, минимизирующего (9), решена. Однако, на практике минимизация и (а) сопряжена с двумя трудностями. Во-первых, вычисление математического ожидания в (6), и, во-вторых, отсутствие во многих экспериментах второго набора данных, т.е. вектора 7*

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

7*

предлагаемый способ формирования вектора у отличается от известных (например, [8]).

Шаг 1. Все измеренные значения

/ = 1... /7. разбиваются на К групп по т измерений в каждой, и в А-ую группу входят измерения:

/к+Ь> /к+2Ь • • ••/\+(т-1)Ь к = \,2.. .К - периодичность выборки при формировании групп измерений. Задавая величину L, можно определить количество измерений в группе и количество групп:

Г п 1

m =

K = L,

где

- целая часть числа.

и вычисляется ве-

личина

L

где

набор

Ik={k,k + L...k + (m-\)L}.

Шаг

3.

После

1 k ( \

(a) = T£u )(a) . (!2)

K

u

k=1

Следует заметить, что вычисленное значение оценки а и будет отличаться от оптимального параметра аор{ (точки минимума функционала (6)) по следующим причинам:

• вычисление функционала и (а) по одной

реализации зашумленного сигнала обуславливает отличие второго слагаемого в (7) от нуля;

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

Все это вызывает теоретическое отличие двух точек минимума: значения а и от оптимального

значения аор1. Поэтому будем рассматривать а и

*

как оценку параметра а ор1, который минимизирует не математическое ожидание ошибки сглажи-

II г II2 вания у — , а ее значение

2

(13)

А(а) = ||у ^п,а для конкретной единственной реализации измеренного сигнала /. По-видимому, значение

a

ор1

для экспериментатора обрабатывающего

(10)

Шаг 2. Из экспериментальных значений исключается ^ая группа, а по оставшимся п — т значениям строится СКС 4—т а )

= 1,2...К, (Н)

индексов

вычисления

иу > (а),А; = 1,2...К , определяется среднее значение:

одну «свою» реализацию зашумленного сигнала является гораздо предпочтительнее, чем значение

а^ - оптимального параметра в среднем для большого числа реализаций.

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

оценок а^ был выполнен обширный вычислительный эксперимент. Остановимся только на некоторых результатах. В качестве тестовых сигналов были приняты два сигнала по спектру: сигнал 1 (низкочастотный) - изображен сплошной кривой на рис. 1 и сигнал 2 (высокочастотный) - показан штриховой линией. Точные значения сигнала, вычисленные в узлах ¿1 < ¿2 <... < ¿п,п = 90 , искажались случайными нормально распределенными погрешностями с относительным уровнем

/"/

. Относительная ошибка фильтра-

5л =

ции

определялась

выражением

Шаг 4. Используя методы одномерной минимизации, вычисляем точку минимума аи функционала (12).

5 f (a) =

n, a

- f

*

к

11

: /; 1 / 1 t vi: i Л -

' s '» \ / -

1---- t

0 0.5 1 1.5

Рис. 1. Тестовые сигналы

В качестве количественной характеристики оптимальности параметра сглаживания а и был введен коэффициент эффективности:

Eu =

s *

n,a

- f

opt

!,aU

- f

(14)

Очевидно, что чем ближе значения коэффициентов к 1, тем меньше проигрыш по точности пара*

метра сглаживания а^ по сравнению с a0pt. Так как этот коэффициент является случайной величиной со значениями в интервале [0,1], то по выборке объемом 200 вычислялись числовые характеристики: среднее значение Ец и минимальное значение min Ец (которое характеризует максимальный проигрыш по точности, который возможен с очень низкой вероятностью - примерно 0.005).

Первоначально рассмотрим случай не коррелированного шума измерения (белый шум). На рис. 2 показаны графики следующих зависимостей, вычисленных для сигнала 1: сплошная - относитель-

ная ошибка 5у (а) при 5^ = 0.10; штрих-точечная - 5у (а) при = 0.03; штриховая - и (а)

при = 0.10; точечная - и (а) при 8^ = 0.03

*

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

вычисления точки минимума функционала и (а).

Особенно этот недостаток заметен для сигнала 2, характеристики которого приведены на рис. 3 (обозначения кривых, как и на рис.2). Из этого рисунка

видно, что при уровне шума 0.03 значения и (а)

(точечная кривая) для а< 10-4 практически постоянны даже в области минимума относительной ошибки сглаживания (штрих-точечная кривая) и такой характер зависимости и (а) может вызвать

большие ошибки вычисления аи, т.е. аи может вычисляться не из области минимума ошибки 5у (а) и существенно отличаться от а0р?.

Рис. 2. Характеристики выбора Рис. 3. Характеристики выбора параметра аи для сигнала 1 параметра аи для сигнала 2

В таблице приведены количественные характеристики Ец, min Ец для сигналов 1,2 (столбцы «Шум не коррелирован»).

Сигнал 1 Сигнал 2

Уровень Шум не Шум Шум не Шум

шума 8Л коррелирован коррелирован коррелирован коррелирован

EU min Ey EU min Ey EU min Ey Eu тт Еу

0.02 0.961 0.717 0.907 0.705 0.899 0.742 Отсутствует однозначное определение величины аи

0.05 0.942 0.628 0.851 0.568 0.942 0.785

0.07 0.926 0.581 0.801 0.491 0.961 0.805

0.10 0.910 0.480 0.736 0.420 0.964 0.835

0.15 0.893 0.438 0.654 0.368 0.966 0.821

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

Рассмотрим случай коррелированного шума измерений. Точные значения сигнала искажались случайными шумами с корреляционной функцией

Я(0) = 1, Я(1) = 0.80, Я(2) = 0.60,

Я(3) = 0.40, Я(4) = 0.20 (сильно коррелированный шум). На рис. 4 приведены графики характеристик (обозначения те же, как и в рис. 2) сглаживания для сигнала 1 при коррелированном шуме измерения, а на рис. 5 - сигнала 2.

Рис. 4. Характеристики выбора параметра Рис. 5. Характеристики выбора параметра аи для сигнала 1 (шум коррелирован) аи для сигнала 2 (шум коррелирован)

Видно, что в случае коррелированного шума существенно изменился характер зависимости

и (а) . Для сигнала 1 точки минимума этого функционала существенно (на 2 порядка) сдвинулись в

сторону меньших значений от величины а,

ор(

а

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

определение а и из-за больших интервалов одинаковых значений функционала и (а) (см. рис. 5) и

в этом случае вычисление этого значения теряет смысл.

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

шума возможно вычисление значения аи суще*

ственно меньшого оптимального значения аор(,

что приводит к недостаточному сглаживанию за-шумленного сигнала, т.е. остается значительный «остаточный» шум. Для ответа на вопрос, как выбрать параметр а в случае коррелированного шума, можно рекомендовать первоначально оценить (например, см. [10]) дисперсию шума, а затем вычислить а на основе соответствующих алгоритмов выбора параметра при известной дисперсии [1,4].

СПИСОК ЛИТЕРАТУРЫ:

1. Воскобойников Ю.Е. Математическая обработка эксперимента в молекулярной газодинамике / Ю.Е. Воскобойников, Н.Г. Преображенский, А.И. Седельников. - Новосибирск : Наука. 1984. -238 с.

2. Сидоров Д.Н. Методы анализа интегральных динамических моделей: теория и приложения / Д.Н. Сидоров. - Иркутск : Изд-во ИГУ, 2013. -293с.

3. Воскобойников Ю.Е. Устойчивые алгоритмы непараметрической идентификации динамических систем / Ю.Е. Воскобойников. - Новосибирск : Изд-во НГАСУ (Сибстрин), 2019. - 156 с.

*

4. Воскобойников Ю.Е. Новый устойчивый алгоритм непараметрической идентификации технических систем / Ю.Е. Воскобойников, В.А. Боева // Современные наукоемкие технологии. - № 5. 2019. - С. 25-29.

5. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач/А.Н. Тихонов, В.Я. Арсенин. -М.: Наука, 1983. - 284 с.

6. Завьялов Ю.С. Методы сплайн-функций / Ю.С. Завьялов, Б.И. Квасов, В.Л. Мирошниченко. -М. : Наука, 1980. - 345 с.

7. Wang Y. Smoothing Splines Methods and Applications. Ser. Monographs on Statistics and Applied Probability v. 121. - A Chapman & Hall book. 2011. - 347 P.

8. Kohavi R. A Study of Cross-Validation and Bootstrap for Accuracy Estimation and Model Selection. — 14th International Joint Conference on Artificial Intelligence, Palais de Congres Montreal, Quebec, Canada. - 1995. - P. 1137-1145.

9. Воскобойников Ю.Е. Устойчивые алгоритмы решения обратных измерительных задач / Ю.Е. Воскобойников. - Новосибирск : Изд-во НГАСУ (Сибстрин), 2007. - 184 с.

10. Воскобойников Ю.Е. Оценивание дисперсий шумов измерений различной статистической природы / Ю.Е. Воскобойников. // Труды НГАСУ. -№ 3. - 2020. - С. 77-90.