Анализ и устранение шумовой компоненты во временных рядах с переменным шагом

Скляр Александр Яковлевич кандидат технических наук доцент, кафедра прикладной математики, Российский технологический университет (МИРЭА) 119602, Россия, г. Москва, пр-т Вернадского, 78 Sklyar Alexander PhD in Technical Science Associate Professor, Department of Applied Mathematics, Russian University of Technology 119602, Russia, g. Moscow, pr-t Vernadskogo, 78 askliar@mail.ru Другие публикации этого автора

1. Введение

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

При анализе данных временного ряда и прогнозировании на его основе возникает множество задач, связанных с выделением трендовой, колебательной и составляющих ^[1,2.3,4]. Для выявления таких рядов используются различные методы в зависимости от характера данных ^[5,6,7,8]. В то же время анализ и обработка зашумленных данных вызывает значительные трудности. Возникает задача устранения, по возможности, такого шума. Для устранения такого шума используются различные методы сглаживания, такие как, методы скользящей средней, экспоненциального сглаживания и др. ^[9], в частности методы, связанные с добавлением белого шума и последующим устранением шумовой компоненты за счет его усреднения ^[10].

2. Оценка шумовой компоненты в исходных данных

Обозначим значения аргумента (временные отрезки) через xi, соответствующие им наблюдаемые значения через y_i, предполагаемую функцию «истинной» зависимости – f(x). И исследуемый ряд представим в виде пар (x_i , f(x_i)). Будем рассматривать общий случай, когда исходные данные задаются с переменным шагом x_i+1=x_i+h_i, где величины h_i, вообще говоря, различны.

Наблюдаемые данные будем представлять в виде , где s_i – шум.

Пусть функцию f(x) имеет производные до 4 порядка включительно, тогда ее значение в точке x+t может быть представлено как

(1)

Рассмотрим значения функции в окрестности точки x₀. Выберем точки x_k=x₀+t_k, где k=0,1,2,3,4; t₀=0 и все t_k различны.

В матричном виде она примет вид.

AB=C

Где ; ;.

Индекс k принимает значения 1, 2, 3, 4.

Ранг матрицы A равен 3, следовательно, существует вектор , где не все λ_k равны 0 такой, что , тогда .

Значения λ_k определяются с точностью до постоянного множителя, в частности, положив λ₄=1, получим допустимый набор из решения системы линейных уравнений

(2)

Определитель матрицы A

Аналогично

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

(3)

В частности, решение этой системы для равноотстоящих узлов ^[11] t₀=0, t₁=-1, t₂=1, t₃=-2, t₄=2 дает значения

Учитывая (2) получаем

(4)

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

Далее учтем, что f(x_m+t_k)=y_m+k-s_m+k и, следовательно

При отсутствии быстрых, то есть с периодами соизмеримыми с шагом ряда, осцилляций f(x) величину можно считать малой, и тогда получаем

(5)

Величина систематических отклонений от 0 будет тем меньше, чем меньше будет , то есть при выборе набора из пяти точек в качестве базовой точки x целесообразно выбирать точку x₃. Перенумеруем точки и введем , тогда (5) примет вид

(6)

Правая часть равенства (6) представляет собой случайную величину. Пусть величины s_m+k – независимые случайные величины с 0 математическим ожиданием и дисперсией σ², тогда математическое ожидание .

Среднеквадратичное значение шума σ2, таким образом, можно оценить исходя из

(7)

3. Выделение в данных функциональной и шумовой компонент

В этих условиях можно определить значения шумовой компоненты s и, следовательно, f(x) исходя из

(8)

Перепишем (7) в матричных обозначениях. Для этого введем матрицу L=(l_ij), где , тогда

(9)

И условие минимума принимает вид

или

(10)

В отличие от случая с равноотстоящими узлами здесь возникает необходимость вычисления элементов матрицы L для каждой строки (при равноотстоящих узлах ненулевые элементы матрицы L одинаковы для всех строк). Система (9) из-за ограничения является нелинейной и нахождение ее решений даже при небольшой размерности встречает значительные вычислительные трудности. В то же время решение системы линейных алгебраических уравнений (СЛАУ) из (10) при заданном значении t не создает особых проблем. В самом деле, матрица является 9 диагональной ленточной матрицей и можно показать, что время решения СЛАУ с такой матрицей линейно по числу уравнений.

Отметим, что матрица является неотрицательно определенной и все ее собственные числа μ_i≥0.

Пусть E_i – собственные вектора матрицы, соответствующие собственным числам μ_i≥0,

Тогда первое из равенств (9) примет вид

В этих условиях представляет собой при t>0 монотонно убывающую функцию от t.

Заметим, что при t→0 минимум функции F(s) в (8) будет достигаться, как следует из (5), при . Последнее условие означает стремление 4-ой производной к 0, то есть исходная функция будет близка к полиному не выше 3 степени.

При t→∞ минимум F(s) будет достигаться, при s_i→0 и, следовательно, функция f(x_i)→y_i.

4. Алгоритм удаления шума и выделения функциональной компоненты в данных

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

1. Вычисляем оценку шума σ² на основе (6). Задаем начальное значение t₀.

2. Решаем СЛАУ в соответствии с (9) и вычисляем значение <s,s>.

3. Если <s,s> > σ², переходим к пункту 4, иначе к пункту 5.

4. Решаем СЛАУ в соответствии с (9) с t=2t₀ и вычисляем значение <s,s>. Если <s,s> >nσ², устанавливаем t₀=t и повторяем пункт 4. В противном случае устанавливаем t₁=t и переходим к пункту 6.

5. Решаем СЛАУ в соответствии с (9) с t=t₀/2 и вычисляем значение <s,s>. Если <s,s> <nσ², устанавливаем t₀=t и повторяем пункт 5. В противном случае устанавливаем t₀=t, t₁=t₀ и переходим к пункту 6.

6. Искомое значение t лежит между t₀ и t₁. Решаем СЛАУ в соответствии с (9) с t= (t₀ + t₁)/2 и вычисляем значение <s,s>. Если |<s,s> -nσ²|<ε, то останавливаем процесс и на основе вычисленных значений s_i рассчитываем значения искомой функции f(x_i)= y_i-s_i. Если нет, то строим новый интервал, устанавливая в зависимости от выполнения неравенства <s,s> >nσ² либо t₀=t, либоt₁=t и переходим к пункту 6.

Отметим, что говоря о шуме и его дисперсии в (6) и (7) предполагается, что он представляет собой одинаково распределенную случайную величину на всем временном ряде. Если это не так, то вместо абсолютной величины шума s_i в (5) шум удобнее представлять в виде s_i=φ(x,y)u_i так, чтобы случайная величина u_i была бы одинаково распределенной на всем временном ряде. Наиболее естественно предполагать при большом разбросе исходных данных, что шум является результатом измерений, которые имеют постоянной относительную погрешность. В этом случае естественно принять s_i=y_iu_i, тогда (6) и (7) примут соответственно вид

(10)

(11)

Тогда дисперсия относительного шума u определяется из (11)

И задача (7) приобретает вид

(12)

Или в матричном виде

Тогда условие минимума принимает вид

или

Отдельно следует отметить, что предлагаемая схема исключает из выделения шума компоненты зависимостей до полиномов 3 степени включительно. Последнее может оказаться обременительным для сильно зашумленных данных. В этом случае для исключения шума удобнее использовать более грубую схему, исключающую из выделения шума компоненты зависимостей до полиномов только 2 степени.

(13)

В частности, решение этой системы для равноотстоящих узлов ^[11] t₀=0, t₁=-h, t₂=h, t₃=-2h, t₄=2h дает значения

5. Результаты численного моделирования

На рисунке 1 представлены результаты обработки зашумленных данных.

Значения по осям x и y – , где rnd представляет собой случайную величину, равномерно распределенную на интервале (-0,1;0,1).

Рисунок 1

На рисунке 2 представлены результаты обработки данных, представляющих значения функции y=e^x на интервале (0;10) с шагом 0,1 округленные до двух значащих цифр.

Среднеквадратичная абсолютная погрешность исходных данных от теоретической кривой в рассматриваемом примере составляет 56, при удалении абсолютного шума - 30, при удалении относительного шума – 19. Таким образом, в случаях, когда данные меняются в широких пределах, а в данном случае отношение максимального значения к минимальному составляет 22000, удаление относительного шума приводит и к лучшему удалению абсолютного шума по сравнению с методом прямого удаления абсолютного шума.

Рисунок 2

6. Выводы

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

· обоснованно определить уровень как абсолютного, так и относительного шума в данных вне зависимости от равномерности шага измерений в исходных данных;

· удалить из данных шумовую компоненту;

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