CC BY
94
УДК 519.246.8+519.254
Вестник СПбГУ. Сер. 1, 2006, вып. 2
Ф. И. Александров
ВЫДЕЛЕНИЕ АДДИТИВНЫХ КОМПОНЕНТ ВРЕМЕННОГО РЯДА ПРИ ПАКЕТНОЙ ОБРАБОТКЕ МЕТОДОМ «ГУСЕНИЦА»-SSA
1. Постановка задачи, история, актуальность
Временным рядом называется последовательность N наблюдений некоторого процесса Ем = (/о, /1,..., /м-1), /п € К, N > 2. Важной задачей является задача выделения аддитивной составляющей ряда, когда требуется найти (или аппроксимировать) Е^ по наблюдаемому ряду Ем = ЕN + Е^ (поэлементно). Существует много вариантов данной задачи: выделение тренда, сезонности, сигнала, — а также подходов для их решения: параметрические и авторегрессионные модели, линейные фильтры, преобразование Фурье, вейвлеты. Среди них можно выделить группу методов, основанных на представлении ряда в некотором базисе. В использованном нами методе «Гусеница»-ЭЭА базис порождается самой структурой ряда, что делает метод гибким и позволяет обрабатывать широкий спектр рядов.
Метод «Гусеница»-88А [1, 2] зародился в 1970-1980-х годах для поиска аттракторов динамической системы, улучшения соотношения сигнал/шум, идентификации сигнала в шуме. В его основе лежит преобразование ряда в матрицу и её сингулярное разложение, используемое в анализе главных компонент. Данный метод позволяет выделять различные аддитивные составляющие ряда, причём выбор нужных для выделения сингулярных векторов делается визуально.
К настоящему времени метод «Гусеница»-ЭЭА хорошо себя зарекомендовал для решения задач анализа и прогноза в различных областях: геофизике, метеорологии, эконометрике, социальных науках (см. [1, 2]). В связи с обработкой значительных объемов данных возникла проблема автоматизации визуальной части метода. В работе [3] были предложены методы автоматического выбора собственных векторов при выделении тренда и периодической составляющей, но не были изучены вопросы выбора их параметров, а также методики их применения. Нами были разработаны подходы к выбору параметров, предложен более устойчивый метод для выделения тренда, а также выработана методика для применения этих методов при обработке серии рядов.
Подобная постановка задачи, в частности, была мотивирована необходимостью применения метода к данными по экспрессии гена, которые представляют собой естественные классы (задаваемые типом гена и возрастом эмбриона), каждый из которых содержит множество незначительно отличающихся рядов (соответствующих конкретному эмбриону).
2. Базовый алгоритм метода «Гусеница»-SSA
Опишем вкратце алгоритм метода «Гусеница»-88А для выделения аддитивной составляющей Ед,-) ряда Ем = Е™ + Е(2) (подробное описание см. в [1]).
На первом этапе алгоритма мы выбираем параметр метода — длину окна Ь, 1 < Ь < N, и строим траекторную матрицу ряда X € Кь хК , К = N — Ь +1, столбцами которой являются вектора X3 = (/—1,..., /]+ь-2)Т,] = 1,...,К. Затем вычисляется матрица ХХТ, её собственные числа {А3-}^=1, ^ ^ ^3+1, собственные вектора {из}^=1, соответ-
© Ф. И. Александров, 2006
ствующие этим собственными числам, где с! = шах{' : > 0}, а также факторные вектора = Хт[/3-/у/А^. Собственные тройки формируют син-
гулярное разложение X = л/^Ч^з^-^■
На втором этапе мы должны идентифицировать группу собственных троек с номерами такую, что она соответствует —(г1'). После этого мы получаем компоненту ряда Р^ (в общем случае приближенно равную с помощью ганкелизации мат-
рицы X1-1) = Х^ед з^з и последующей развертки её в ряд обратно тому, как мы
делали это при построении траекторной матрицы ряда.
Также можно прогнозировать выделенную компоненту , так как «Гусеница»-SSA позволяет найти линейную рекуррентную формулу, которой подчиняется компонента.
Вопрос выбора Ь подробно описан в [1, 4]. При заданном Ь формирование результата управляется выбором собственных троек Д. Количество необходимых троек зависит от вида —Например, гармонике (синусу с произвольным фазовым сдвигом) с периодом Т > 2 соответствуют две тройки, полиному степени т соответствуют т + 1 троек, а экспоненте — одна (см. теорию рядов конечного ранга в [1, 4]). Общее правило по выбору троек состоит в том, что при выполнении некоторых ограничений ряду определённого вида соответствуют тройки с собственными векторами, последовательность элементов которых имеет тот же вид. Например, гармонике соответствуют две тройки с векторами, последовательности элементов которых задаются формулой синуса/косинуса с той же частотой (причём собственные числа этих троек равны), полиному степени т соответствуют вектора, последовательности элементов которых задаются полиномами порядка не выше т.
Таким образом, тренду, определённому как медленно меняющаяся составляющая ряда, соответствуют «медленно меняющиеся» собственные вектора, а каждой гармонике, составляющей периодическую составляющую ряда, соответствуют два собственных вектора вида синус-косинус с одинаковой частотой. Пользуясь этими правилами, можно формализовать поиск необходимых собственных троек, перейдя на язык разложения Фурье.
3. Критерии для идентификации
собственных троек тренда и гармоники
Будем рассматривать по очереди все собственные вектора полученного разложения ряда —N. Будем считать периодограмму последовательности элементов каждого собственного вектора V = (и\, и^,.и^)Т со следующей нормировкой:
т . 2со2, к = 0,
2с^/22, если Ь — чётное и к = Ь/2,
где вк и в к —это коэффициенты разложения Фурье для последовательности и\, и2,.. -,иь при косинусе и синусе с частотой к/Ь. Значение П^ (к/Ь) отражает вклад гармоники с частотой к/Ь в разложение Фурье последовательности элементов V.
Если собственный вектор некоторой тройки таков, что гармоники с низкими частотами имеют большой вклад в разложение Фурье последовательности его элементов (это соответствует тому, что она медленно изменяется), то отнесем эту тройку к трендовым.
Это можно проверить следующим образом:
(к/Ь) ^
с\и) = ^-т-гЫи/тл ^ 0
2^о<к/ь<о.б пи(к/Ь)
для заданного Со € (0,1]. Задавая и>о, мы определяем область низких частот [0,о>о].
Для выделения гармоники (с Т > 2) будем искать два вектора из троек с номерами 2,2 + 1, периодограммы которых сосредоточены в одной и той же частоте. Будем считать, что это так, если
Pjj+1 = 0.5 o<max (п£. (k/L) + п£.+1 (k/L)) > ро
для заданного ро. Для неискаженной гармоники p¿,¿+i должно быть равно 1, также было теоретически получено асимптотическое значение p¿,¿+i для экспоненциально-модулированной гармоники (см. [5]).
Проведенные с помощью моделирования исследования показали, что при оптимальных пороговых значениях Со, ро данные критерии позволяют с хорошим качеством выделить гармонику из суммы гармоники и белого нормального шума, и тренд из за-шумленного полинома (экспоненты).
Применение методов к серии реальных рядов ставит перед нами два вопроса. Во-первых, работают ли методы на рядах данной серии, т. е. достигается ли достаточно хорошее качество и, во-вторых, какие пороговые значения надо выбирать.
Оценка качества методов для обработки серии рядов. Покажем, как можно оценить качество на примере выделения тренда. Возьмём тестовое подмножество, для каждого ряда Fn из этого подмножества с помощью визуального исследования выделим искомую составляющую F^. Для этого ряда найдём такое пороговое значение C0pt, которое обеспечивает наилучший результат —^(С^) по сравнению с «визуальным», и вычислим норму разности между ними: HF^ — (C^pt) Уг2. Среднее значение её по тестовому подмножеству считаем оценкой качества в среднем по всей серии. Если оно достаточно мало, то считаем, что метод на рядах из серии работает. Размер тестового множества зависит от того, насколько похожие (по виду, а, следовательно, по структуре составляющих их собственных троек) ряды составляют всю серию.
Выбор оптимальных пороговых значений. Выбор оптимальных ро, Со проводится по-разному, так как в первом случае известен параметрический вид искомой составляющей, что упрощает задачу. Мы можем посчитать ро, исходя из данных об уровне модуляции гармоник, характере шума, отношении сигнал/шум. Зная это, мы можем с помощью моделирования получить искомое ро (см. [5]).
Для тренда же, ввиду неизвестной нам параметрической модели, мы должны подбирать Со, исходя из вида данных, в идеале следующим образом: Со = argminc0 II F^ — —(1)(Со))||г2. Но F(1) неизвестно, и нами была предложена другая мера качества, основанная на том факте, что при уменьшении Со в F(1) попадают компоненты, у которых П^ всё меньше сосредоточена в [0,о>о]:
V(F Ж )]_C(FN)-C(FN-F^)
Мера R(Fn , Fn\Co)) при изменении C0 ведёт себя подобно \\Fff — F^ (Co))||i2 и не зависит от Fn^- Таким образом, для любого ряда из серии можно найти Co, при котором результат метода будет наиболее близок к визуально выделенному результату:
Co = argminj R(Fn ,F(1)(Co))}.
Разработанная методика автоматической идентификации была применена к модельным рядам, имитирующим реальные зашумленные ряды. Модели тренда и шума были взяты на основе результатов исследования уровня экспрессии гена Kruppel мушки дрозофилы. Результаты оказались лишь чуть хуже результатов визуальной идентификации, которая в данной постановке задачи считается наилучшей.
Summary
F. I. Alexandrov. Batch extraction of additive components of time series by means of the «Caterpillar»-SSA method.
Criteria for eigentriples identification in «Caterpillar»-SSA are offered. They allow one to automate the extraction of trend and periodical components of time series. Approach to choice of parameters of criteria is given and technique for batch processing based on these criteria is described.
Литература
1. Golyandina N., Nekrutkin V., Zhigljavsky A. Analysis of Time Series Structure: SSA and Related Techniques. Chapman&Hall/CRC, 2001. 305 p.
2. Ghil M, Allen R.M., Dettinger M.D., Idee K, Kondrashov D, Mann M.E., Robertson A., Saunders A., Tian Y., Varadi F., Yiou P. Advanced spectral methods for climatic time series // Rev. Geophys. 2002. Vol. 40, N 1. P. 1-41.
3. Vautard R., Yiou P., Ghil M. Singular-Spectrum Analysis: A toolkit for short, noisy chaotic signals // Physica D. 1992. Vol. 58. P. 95-126.
4. Голяндина Н. Э. Метод «Гусенип^-SSA: Анализ временных рядов. СПб., 2004. 76 с.
5. Alexandrov Th., Golyandina N. Automatic extraction and forecast of time series cyclic components within the framework of SSA // Proc. of the 5'th Worksh. on Simul. St.Petersburg, 2005. P. 45-50.
Статья поступила в редакцию 20 мая 2006 г.
