CC BY
64
_Доклады БГУИР_
2010 № 2 (48)
ИНФОРМАТИКА
УДК 001:002
ИССЛЕДОВАНИЕ SSA-МЕТОДА НА ОСНОВЕ КОМПЛЕКСНОГО ПРИМЕНЕНИЯ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ
Е.А. ЯКИМОВ
Белорусско-Российский университет пр. Мира, 43, Могилев, 212005, Беларусь
Поступила в редакцию 7 декабря 2009
Предложена методика исследования SSA-метода с комплексным применением информационных технологий. Экспериментальные исследования детерминированного временного ряда проведены с использованием табличного процессора MS Excel, математического пакета MathCAD и пакета статистической обработки данных Statistica. Определены количественные характеристики погрешности восстановления временного ряда.
Ключевые слова: SSA-метод, временная последовательность данных, MS Excel, MathCAD, Statistica.
Введение
В процессе интеллектуального анализа данных (ИАД) центральное место занимает автоматическое порождение характеризующих анализируемые данные моделей, правил и/или функциональных зависимостей. В целом процесс извлечения знаний в ИАД условно делят на следующие этапы [1], которые в совокупности предложено использовать на этапе эксплуатации имитационной модели (ИМ) сложного объекта [2].
Шаг 1. Отбор данных: анализ задач пользователя, выбор целевого множества данных, определение переменных.
Шаг 2. Предобработка данных: устранение зашумленности, обработка пропущенных значений, итоговые показатели по группам данных.
Шаг 3. Редукция и проекция данных: ищутся полезные особенности данных для решения поставленных задач, сокращается пространство переменных.
Шаг 4. Поиск закономерностей: выбор метода поиска закономерностей с учетом объема и типа данных, их зашумленности и осуществление поиска закономерностей.
Шаг 5: Оценка и интерпретация найденных закономерностей: оценка и упорядочение закономерностей по их релевантности, проверка согласованности предыдущих и вновь найденных знаний. Возможно возвращение к любому шагу от 1 до 4 для дальнейших итераций.
Шаг 6. Использование найденных знаний: прямое использование, передача заинтересованным лицам, включение в интеллектуальные системы, основанные на знаниях.
Для разработки технологии извлечения знаний из временных последовательностей данных исследован сингулярный спектральный метод (SSA-метод), включающий этапы вложения, сингулярного разложения, группировки, диагонального усреднения [3].
Содержательное описание SSA-метода
SSA-метод используется для анализа временных рядов и может быть использован на каждом из этапов эксплуатации ИМ. Он позволяет выделить ряды, которые описывают тренд
первоначального ряда, гармонические колебания и те составляющие рядов, которые относят к "шуму". При этом метод не требует стационарности ряда, знания модели тренда, а также сведений о наличии в ряде периодических составляющих и их периодах. Также с помощью данного метода можно определить модель тренда и использовать это знание для дальнейшей обработки ряда уже с известной моделью тренда, что важно, например, при автоматизации определения длительности переходного процесса в ИМ.
Математической основой SSA-метода является сингулярное разложение [4]. Для успешного применения SSA-метода следует последовательно пройти несколько шагов.
Вложение. На этом шаге выбирается ширина окна L, от выбора которой зависят результаты применения SSA-метода. Из-за того, что нет общих рекомендаций по выбору ширины окна, параметр L зависит от решаемой задачи и предварительной информации, известной о ряде. Например, для выделения тренда рекомендуется выбирать ширину окна не слишком большой. С другой стороны, для выделения гармонических колебаний рекомендуется большая ширина окна. После выбора ширины окна в соответствии с L строится траекторная матрица А ряда, которая будет являться по условию ее построения ганкелевой [5].
Сингулярное разложение. Для матрицы S=A-AT находятся собственные числа X и орто-нормированные собственные вектора U. Упорядоченные по убыванию собственные числа, которые являются большими нуля, часто называются сингулярными числами, а соответствующие им собственные вектора — левыми сингулярными векторами U. После этого вычисляются вектора V, которые называются правыми сингулярными векторами, и находятся элементарные матрицы, на сумму которых раскладывается первоначальная траекторная матрица.
Группировка. На этапе группировки элементарные матрицы группируются по принципу принадлежности к тренду, гармоническим колебаниям или к шуму. Этот этап является наиболее сложным при применении SSA-метода. Для нахождения тренда на диаграммах собственных векторов (по оси абсцисс откладывается порядковый номер координаты собственного вектора, по оси ординат откладывается значение координаты собственного вектора) выделяют медленно меняющиеся вектора. Сумма элементарных матриц, соответствующих этим векторам, будет являться траекторной матрицей тренда ряда. После этого восстанавливают гармонические колебания ряда. Для отделения шума можно воспользоваться несколькими замечаниями: нерегулярное поведение сингулярных векторов может говорить о принадлежности их к набору, порожденному шумовой компонентой; также об этом может свидетельствовать медленное, практически без скачков, убывание собственных значений с некоторого номера.
Диагональное усреднение. Если полученные сгруппированные матрицы являются ганке-левыми, то они являются траекторными матрицами некоторого ряда, который может быть легко по ним восстановлен. Однако обычно сгруппированные матрицы редко получаются ганкеле-выми, поэтому для восстановления ряда прибегают к диагональному усреднению. В соответствии с этим этапом каждый член восстановленного ряда будет являться средним арифметическим соответствующей ему побочной диагонали траекторной матрицы.
В результате проделанных шагов получается несколько рядов, один из которых является рядом, описывающим тренд первоначального ряда, другой описывает гармонические колебания, а третий — шумовые составляющие.
Методика исследования SSA-метода на основе информационных технологий
Для исследования SSA-метода применяется комплекс информационных технологий, представленный табличным процессором MS Excel, математическим пакетом MathCAD и пакетом статистической обработки данных Statistica.
Этап вложения. Для экспериментальных исследований исходный ряд G=(g0, gi, ..., gN-i) может быть задан по известным функциям либо определен результатами функционирования исследуемой системы. По известным функциям ряд формируется на рабочем листе MS Excel и затем в пакете MathCAD формируется матрица А, которая по правилам построения является ганкелевой [5]. Процедура вложения является преобразованием исходного одномерного ряда G=(go, gi, ..., gN-i) в последовательность L-мерных векторов, число которых равно K=N-L+i:
A=(gi. J , g+. J, 1J i J K . (i)
Эти вектора образуют траекторную матрицу Л=[А^ ...; А^] ряда О, в которой aij=gi+j_2, т. е. матрица А имеет одинаковые элементы на диагонали /+/—сопз!.
Этап сингулярного разложения. Обозначим Б=АА Матрица А А симметричная и
неотрицательно определенная, а значит ее собственные числа \к ^ вещественны и неотрицательны. Представленные в виде /. | >... >/.,>() собственные числа называют сингулярными значениями матрицы А. Пусть III, ..., иь — соответствующие им ортонормированные собственные вектора. Будем называть /;=тах {к | /./>()} порядком сингулярного разложения. Обозначим
1
V, = -= А и,, к = 1,...,р . (2)
V \
Тогда сингулярным разложением матрицы Л называется ее представление в виде суммы элементарных матриц
А = А1+А2 + ... + Ар, Ак=^Тк Мк\тк
Каждая из матриц Ak имеет ранг, равный единице. Поэтому их можно назвать элементарными матрицами. Вектор U^ называют к-м левым сингулярным вектором или просто к-м собственным вектором, вектор V/. — правым сингулярным вектором. Набор <V/Wl.. U^, V/> будем называть к-й собственной тройкой.
Собственные числа -\.к ^ матрицы А в пакете Math CAD представлены вектором d.
Вектор d сингулярных значений в MathCAD определяется с использованием функции svds() [6]:
d := svds(A). (4)
Диагональная матрица ds сингулярных значений матрицы А в пакете MathCAD определяется с использованием функции diag():
ds := diag(d). (5)
Объединенная матрица AS с левыми и правыми сингулярными векторами определяется с использованием функции svd ():
AS := svd(A). (6)
Для разделения левых и правых сингулярных векторов из матрицы AS используется функция submatrix() [6].
Этап группировки. Вид левых и правых сингулярных векторов, трактуемых в SSA как временные ряды, является очень важным для следующего шага метода — группировки [3]. При этом для одномерного SSA левые и правые сингулярные векторы обладают определенной симметрией, так как в этих случаях сингулярные разложения траекторных матриц с длиной окна L и K=N-L+1 эквивалентны.
Процедура группировки формально одинакова для всех разновидностей SSA. На основе разложения (3) процедура группировки делит все множество индексов р на m непересе-
кающихся подмножеств I1,..., Im .
Пусть I={/'i, ..., ip} Тогда результирующая матрица AI, соответствующая группе I, определяется как A = A +K +A . Такие матрицы вычисляются для I=I1, ..., Im, тем самым разложение (3) может быть записано в сгруппированном виде:
A=A +K +A . (7)
1 m
Процедура выбора множеств I=I1, ..., Im и называется группировкой собственных троек. Для определения I=I1, ..., Im в MS Excel используется лепестковая диаграмма, которая является аналогом графика в полярной системе координат, отображая распределение значений относи-
тельно начала координат. По особенностям представления сингулярных векторов на лепестковой диаграмме принимается решение о принадлежности их одной группе.
Этап диагонального усреднения. На последнем шаге базового алгоритма каждая матрица сгруппированного разложения переводится в новый ряд длины N. Для произвольной матрицы X процедуру приведения ее к ганкелевому виду и последующему преобразованию в ряд (обозначим его как С ) выразим следующим образом. Пусть X — матрица размера ЬхК с элементами хгу, 1 </</,, 1 </<К. Положим Ь =шт(Ь,К), К* =шах(£, К) и К=^Ь+1. Пусть ^у=Ху, если Ь<К и г,=хр в остальных случаях. Тогда диагональное усреднение переводит матрицу X в ряд
по формуле
м k+1 lkе zhk_ + 2 | 0 J k J L* - 1;
j= i
L
gl = н tt е zj
>k H L е j,k- j+ 2 j= 1
N- K* +1
IL* - 1J k J K*;
(8)
nr е Zjik. j+ 2 IK j kj N.
j= k- K + 2
Это выражение соответствует усреднению элементов матрицы вдоль побочных диагоналей i+j=k+2: выбор k=0 дает gB0=x1b для k=1 получаем gB1=(x12+x2i)/2 и т. д. Применив диагональное усреднение к матрицам, полученным на этапе группировки, приходим к разложению исходного ряда в сумму m рядов.
Процедуру диагонального усреднения просто и наглядно предложено выполнить в MS Excel. Для этого матрица, подлежащая диагонализации, размещается на рабочем листе. Затем блок матрицы, следующий за первой строкой, сдвигается вправо на одну позицию. В сдвинутом блоке также определяется блок, следующий за первой строкой, который сдвигается вправо на одну позицию. Процедура повторяется до тех пор, пока в очередном блоке не останется ни одной строки. Восстановленный ряд G определяется аналогично формуле (1) с использованием функции СРЗНАЧ() в MS Excel. Затем исследуется в пакете Statistica.
Экспериментальная часть
Для проведения исследований выбрана функция F(x)=y1(x)+y2(x), где y1(x) — линейная функция, y1(x)=0,1x+2; y2(x) — гармоническая функция, y2(x)=sin(x). Переменная x принимает значения от 0 до 42 с шагом, равным единице. Таким образом, длина N ряда G=(g0, g1, ..., gN-1) равна 43. При этом длина окна L=22, число L-мерных векторов K=N-L+1=22.
Отдельно на рабочем листе MS Excel 2003 рассчитаны значения функций y1(x), y2(x) и F(x). При этом ряду Gyi линейной функции, или тренда, соответствуют значения Gy\ = {gyi0 -Яг!,--- Syi42} = {2,°; 2,1;...; 6,2} . Ряд Gy2 гармонической составляющей —
Gy2={gy2o,gy2l,...,gy2J = {0-,0,m-,...-,-0,9n} HG = {go,gi,...,g42} = {2,0;2,94;...;5,28}.
Элементы ряда G копируются в траекторную матрицу A на рабочем листе MathCAD 14.0. Сформированная матрица A является ганкелевой размером 22^22. Следует отметить, что размеры матрицы при ручном вводе ограничены возможностями MathCAD и не должны превышать 600 позиций.
Сингулярные числа при сингулярном разложении траекторной матрицы А с применением функции svds() по формуле (4) принимают следующие значения: /ч =99.49: л2= 10.99: ^з=10,80; /.4=2.10: /.^=0.024. Остальные значения /.,,. ..., Х22 представляют убывающую последовательность: Х6, ... Д22=0,022; 0,018; ..., 3,06-10~3; 6,602-10^.
Сумма всех сингулярных чисел равна 116,577, причем сумма собственных чисел ^+^+...+^5=116,413.
Сингулярное разложение в пакете MathCAD 14.0 реализуется с помощью функции svd() по формуле (6), формируется обобщенная матрица AS, включающая левые и правые сингулярные вектора. Левые и правые сингулярные вектора разделяются по формулам U:=submatrix(AS, 0,21,0,21) и F:=submatrix(AS, 22,43,0,21), соответственно.
Диагональная матрица ds сингулярных чисел определена по формуле (5). Для контроля правильности выполненных вычислений определяется равенство
U -ds -VT = А (9)
Для последующих расчетов определяются скалярные значения собственных чисел по формулам с использованием функции submatrix(). Например, для /., (обозначение dsO в пакете MathCAD)
ds0 := submatrix (ds, 0, 0, 0, 0), ds0 = 92,49, (10)
для X2 (обозначение ds 1 в пакете MathCAD) ds\ := submatrix(fife, 1,1,1,1), сЫ = 10,99 ит д
Левые сингулярные вектора определяются из матрицы U аналогично с использованием функции submatrix():
U0 := submatrix (U, 0, 21, 0, 0), (12)
и ••= submatrix(U, 0, 21,1,1) и т. д. (13)
Правые сингулярные вектора определяются в пакете MathCAD по формуле 1
V •=-— ЧА 4U,., i = 0,1,..., 4. (14)
г yjds0
Элементарные матрицы X определены по формуле X := /dsiЧU 4VT), i = 0,1,..., 4 . (15)
Для определения групп элементарных матриц X] построены лепестковые диаграммы для векторов V в табличном процессоре MS Excel 2003 (рис. 1 и 2). Вектора Vb V2 и V4 имеют вид, аналогичный изображению на рис. 2. Векторам У0 и У3 на диаграммах соответствуют расходящиеся кривые (рис. 1).
Рис. 1. Лепестковая диаграмма вектора У2 Рис. 2. Лепестковая диаграмма вектора У3
Тогда группа XI элементарных матриц соответствует тренду ряда: Х1:=Х0+Х3. Аналогично группа Х2:=Х1+Х2+Х4 соответствует гармонической составляющей ряда.
На этапе диагонализации для выполнения расчетов по формуле (8) используется специальная процедура в MS Excel, которая иллюстрируется на примере матрицы X1 с левым верхним блоком
ж2,141 2,242 2,276 2,3034
i2,242 2,345 2,377 2,4024
з2,276 2,377 2,407 2,431ч
Й2,303 2,402 2,432 2,4531
На рис. 3 представлен фрагмент рабочего листа MS Excel с решением задачи диагонализации по предложенной процедуре.
Сдвиговая процедура блока матрицы X1 при диагонализации 2,141 2,242 2,276 2,303 x115 Х116 x117
2,242 2,345 2,377 2,402 x125 x126
2,276 2,377 2,407 2,431 X135
2,303 2,402 2,431 2,453
Восстановленные значения тренда, g ^ 2,141 2,242 2,299 2,34
Значения тренда, g 2 2,1 2,2 2,3
Отклонения g - g в -0,141 -0,142 -0,099 -0,04 -0,0037 0,069 0,147
Рис. 3. Процедура диагонализации матрицы X1 тренда в MS Excel Восстановленные значения тренда g еу1 определяются по столбцам с применением функции СРЗНАЧ(). Полученная таким образом матрица используется в пакете Statistica 6.0 для оценки характера распределения ошибки восстановления ряда.
Результаты и их обсуждение
В ходе проведенных исследований временного детерминированного ряда G, образованного функцией вида F(x) = ax + 6sinx + c = 0,1x + sinx + 2, построена ганкелева матрица А размером 22^22. При сингулярном спектральном разложении матрицы A основными для восстановления ряда явились элементарные матрицы Х0 и Х3 для линейной составляющей (собственные числа /.| =99.49 и Л4=2.10. которым в пакете MathCAD соответствуют обозначения dsi) и ds3) и матрицы Хь Х2. Х4 для гармонической составляющей (собственные числа /-2= 10.99: ^з=10,80; ^5=0,024, которым в пакете MathCAD соответствуют ds 1. c/s2 и ds4).
Абсолютная ошибка восстановления ряда определена величиной отклонения
D= max\gt - g' \. (17)
Относительная ошибка восстановления ряда определяется по формуле
I в I
d= max\g - ge\ 100% . (18)
max \ g, \
Абсолютная ошибка восстановления для тренда составляет 0,142, для гармонической составляющей — 0,143, для исходного ряда — 0,0055. Относительная ошибка восстановления тренда составляет 2,29%, для гармонической составляющей — 14,3 %, для исходного ряда — 0,08 %.
Обнаружены интересные факты: абсолютная погрешность восстановления тренда и гармонической составляющей для начальных и конечных значений ряда имеют наибольшие значения (рис. 4), распределения суммы квадратов отклонений значений побочных диагоналей матрицX1 и X2 имеют характерный вид (рис. 5).
Рис. 4. Распределение абсолютной ошибки Рис. 5. Распределение суммы квадратов
восстановления тренда отклонений значений побочных диагоналей XI
Для тренда наибольшая ошибка принимает отрицательные значения, для гармонической составляющей — положительные значения. Выявленные особенности оказывают существенное влияние на прогнозирующие свойства 88Л-метода.
Заключение
Предложена методика исследования 88Л-метода, позволяющая проводить его изучение без разработки специализированного программного приложения. Получены оценки погрешности восстановления временной последовательности данных для детерминированной функции Г(х) = ах + ^пх + с.
Работа выполнена при финансовой поддержке гранта Ф09М-171 БРФФИ.
RESEARCH SSA-METHOD ON THE BASIS OF COMPLEX USING INFORMATION TECHNOLOGY
E.A. YAKIMAU Abstract
The technique of research of a SSA-method with complex application of information technologies is offered. Experimental researches determined time of some are lead with use of tabulared processor MS Excel, mathematical package MathCAD and a package of statistical data processing Statis-tica. Quantitative characteristics of an error of restoration time are certain of some.
Литература
1. Таран Т.А., Зубов Д.А. Искусственный интеллект. Теория и приложения: учеб. пособие. Луганск, 2006.
2. Якимов Е.А., Максимей И.В. // Новые материалы, оборудование и технологии в промышленности: материалы междунар. науч.-техн. конф. молод. ученых. Могилев, 19-20 ноября 2009 г. Могилев, 2009. С. 121.
3. Голяндина Н.Э. Метод 'Тусеница''-SSA: анализ временных рядов: Учебное пособие. СПб., 2004.
4. Голуб Дж., Ван Лоун Ч. Матричные вычисления: пер. с англ. М., 1999.
5. Гантмахер Ф.Р. Теория матриц. М., 1966.
6. Ивановский Р.И. Компьютерные технологии в науке и образовании. Практика применения систем MathCAD Pro: учеб. пособие / Р.И. Ивановский. М., 2003.
