CC BY
53
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
Изв. вузов «ПНД», т. 14, № 6, 2006 УДК 530.182
СИСТЕМА ИТЕРАТИВНЫХ ФУНКЦИЙ И МАРКОВСКИЙ ПРОГНОЗ ВРЕМЕННЫХ РЯДОВ
Н.Г. Макаренко, Л.М. Каримова, С.А. Мухамеджанова, И.С. Князева
В статье изложена методика вероятностного прогноза временных рядов на основе системы случайных итеративных функций из теории фракталов. Итерации приводят к аттрактору (фракталу) в пространстве компактов. Аттрактор является носителем инвариантной вероятностной меры (мультифрактала) в пространстве борелевых мер. Обратная задача состоит в нахождении системы итеративных функций и их вероятностей по оценкам эмпирической меры. Такие оценки можно получить из временного ряда, используя методы символической динамики. Кроме необходимых математических сведений, мы приводим пример практического предсказания пороговых значений геомагнитных возмущений.
Решающим для всякой идеи является не то, как она осуществляется, но что, по существу, в ней содержится.
Стефан Цвейг
Введение
Теоретические основы предсказания на основе системы итеративных функций (iterated function system, IFS) опираются на элементы теории меры и некоторые разделы фрактальной геометрии, выходящие за рамки пособий и монографий по теории фракталов, популярных среди большинства физиков. Заполнить этот пробел мы попытались в настоящей статье1. Наша главная цель - резюмировать математические сведения о фундаментальной теореме IFS и вещах вокруг нее на уровне технической строгости. Поэтому большая часть статьи имеет методический характер. В последнем разделе мы иллюстрирует теорию нашими оригинальными результатами по предсказанию магнитных бурь.
Современные методы детерминированного прогноза временных рядов основаны на реконструкции топологической модели динамической системы в евклидовом пространстве Rm. Для этого, предполагается, что система описывается пото-
1Основой послужила лекция первого автора на XIII школе-конференции «Нелинейные волны-2006», Нижний Новгород, 1-7 марта 2006.
ком ф*(у) : М — М на п-мерном компактном2 многообразии М. Мы наблюдаем проекцию траектории ф*(уо) как гладкую морсову функцию3 фазовых точек: хг = Н (у (гх)) : М — К и получаем временной ряд {хг} , хг = х (гх), г = 0,1,... Тогда, отображение Г : М — Кт
Г (у) = (Н (у) ,Н (ф(у)),...., Н(фтт (у)),
с точностью до предположения о типичности, является топологическим вложением Г (М) с Кт, если т = 2п + 1 [1, 2]. В общем случае, когда фазовое пространство М имеет произвольную размерность, для вложения требуется существование п-мерного аттрактора [3]. Состояние полученной модели задается вектором хг = {хи,х(г+1)х, ...,Х(г+т-1)х}, составленным из отсчетов ряда в пространстве вложения Кт, а динамика описывается уравнением Хг+1 = Ф (хг), где Ф = ГофОГ-1. Образ вложения Г (М) является гладким подмногообразием в Кт и диффеомор-фен М. Следовательно, полученная копия наследует все топологические характеристики оригинала [2-4]. Предиктором служит единственная нетривиальная компонента хг+т = Ф (хг) векторного уравнения, а предсказание заключается в поиске подходящей аппроксимации Ф [5, 6], относительно которой известно лишь, что это непрерывная и, возможно, дифференцируемая функция т переменных, определенная на конечной «обучающей» выборке {хг+т, Ф (хг)}. При этих условиях задача аппроксимации не является корректной [6, 7] и решается только на уровне математически правдоподобных утверждений с помощью локальных [8-10] или глобальных методов [6, 11]. Нестационарность временного ряда и присутствие шумов существенно влияют на ошибку модели, выбор которой является сложной проблемой [5, 11-13].
Вероятностное предсказание используют в более общих ситуациях, когда не выполняется кредо идеального экспериментатора [3], необходимое для корректной реализации описанного подхода. Для предиктора чаще всего используют стохастические модели, основанные на марковских процессах [10, 14].
Рассмотрим простую авторегрессионную модель: хп+1 = Тхп = (ахп + ^п), где случайная переменная подчиняется гауссову распределению N о) с нулевым средним и дисперсией о2. Эволюция точки хп — хп+1 полностью описывается условной переходной вероятностью [10]
Р (хп+1 \хп ) = Р (хп+1,хп)/р (хп),
где Р (хп+1,хп) = /р (хп) N (I, о)б (хп+1 - ахп - Ы ^ = N (хп+1 - ахп, о). Если следовать принципу марковских цепей4 - будущее зависит от прошлого только через настоящее, р (хп+1 \хп ,хп-1,.....,х0) = р (хп+1 \хп). Этот пример иллюстрирует идею случайной динамической системы [15, 16]. Последняя определяется набором непрерывных преобразований Т = {Тк} ,Тг : X — X, к = 1,М, заданных на компактном подмножестве X е Кп так, что каждое Тк выбирается с помощью некоторого случайного процесса. Пусть о^ о2, о3,... е {1, 2,.., N} - последовательность символов, снабженных вероятностями рп = Р (оп = п), а хо £ X -
2Компактом называют ограниченное и замкнутое множество.
3Гладкая морсова функция либо не имеет критических точек, либо все они изолированные и невырожденные.
4Цепью называют марковский процесс с дискретным временем.
начальная точка. Тогда динамика случайной системы описывается уравнениями вида: хп+1 = Т0п+1 хп. Мы получим марковский процесс, если каждый текущий выбор отображения зависит от предыдущего, то есть рп = Р (оп = п) = рп,0п_1. Состояние случайной системы (X, Т, Р) можно идентифицировать лишь в некотором грубом приближении. Пусть {Л1,Ап} ,Х С у Лг - разбиение пространства X на конечное число клеток Лг. Игнорируя эволюцию внутри клетки, будем считать, что система имеет состояние5 г, если ее орбита находится в клетке Аг. Марковская динамика полностью описывается (п х п) - матрицей вероятностей Р^ = РгоЬЦ {Тх е Л2 \х е Лг}, где Р^ - вероятность перехода г — ]. Очевидно, что вероятности определятся долей тех точек Лг, которые попали в Л2 под действием Т
Ргз = Ц (Лг П Т-1Л3)/ц (Лг).
Здесь ц(В) инвариантная мера, которую можно определить для почти каждой случайной орбиты Х1,Х2,...., XI и множества В С X, как предел величины #(В П х1,х2,...., XI)/1, где числитель - число точек орбиты, которые попали в множество В, а знаменатель I - длина орбиты [16].
Вероятностное предсказание с помощью марковской модели сводится к двум задачам.
Первая задача заключается в построении грубого разбиения состояний системы на основе наблюдаемого временного ряда для эмпирической оценки инвариантной меры. Это удобно сделать, используя методы символической динамики [17]. Идея заключается в преобразовании ряда в последовательность символов из конечного алфавита 2 = {«1, ..,8п}, кодирующего крупнозернистое разбиение X. Проще всего преобразовать временной ряд {хг} в бинарный «текст», состоящий из слов Б = з1,82,...8к, где зг е {0,1}; обобщение на п > 2 тривиально. Выберем некоторое пороговое значение ординаты хг = Н и трансформируем все отсчеты ряда в символы по правилу: 8г = 0, если хг < Н и 8г = 1, если хг > Н. Разделим единичный отрезок [0,1) на 2К интервалов и каждому из 2К возможных слов длины К поставим в соответствие полуоткрытый сегмент [х[,хг), определенный двоичным представлением слова Б [17, 18]
К / /
х1 (Б) = Е- ,8г/2г; хг (Б) = х1(Б) + 1/2К.
Теперь подсчитаем частоту встречаемости различных слов Б на единичном отрезке. Если полученная гистограмма (1 (х) стационарна6, то ее можно считать оценкой инвариантной меры ц (х). Для многих процессов полученная таким способом мера имеет мультифрактальные свойства [18-20].
Второй задачей марковского предсказания является построение теоретической модели, в которой мера генерируется случайной динамической системой, реализованной набором сжимающих линейных отображений {Т^}, снабженных вероятностями. Такие отображения хорошо известны в геометрии фракталов. Они называются гиперболическими1 IFS [21, 22] и были введены в 1985 году М. Барнсли
5Мы отождествляем состояние г с г-м моментом времени.
6То есть гистограммы для двух фрагментов ряда практически совпадают.
7 Набор линейных отображений подобия {fi} = 1,N,fi : (Х,й) ^ (Х,й) таких, что й (^ (х) , ^ (у)) = rid (х, у) называют гиперболическим, если Уг, \г-1\ < 1.
и С. Демко [23], но, фактически, уже использовались в 1981 году Дж. Хатчинсоном [24]. IFS рассматривают фракталы как предельные образы - аттракторы коллективной непрерывной динамики линейных отображений [21-24]. Таким образом, исключается столь непривычное в континуальной математике «хирургическое» удаление фрагментов геометрического объекта для построения некоторых фракталов. Например, канторово множество K является пределом бесконечного числа итераций объединения двух отображений T1x = (1/3)x и T2x = (1/3)x + 2/3, приводящих к самоподобному и фрактальному множеству K = T1 (K) U T2 (K). Если снабдить IFS вероятностями, то ее динамика продуцирует мультифрактальную меру на аттракторе [21], которая является единственной и инвариантной относительно действия марковского оператора [21, 23, 25]. Случайные IFS нашли интересные применения для генерации и сжатия изображений [21, 25, 26], моделирования геномов [20, 27] и анализа текстов в лингвистике [28, 29]. Если известен аттрактор и мера на нем, то обратная задача в теории IFS заключается в поиске коэффициентов отображений и их сопутствующих вероятностей [30]. Предположим, что мы имеем дело с бинарной символической последовательностью слов фиксированной длины, составленных из символов {0,1}. Если относительная частота встречаемости различных слов стационарна и статистически самоподобна, полученную гистограмму можно рассматривать как эмпирическую оценку мультифрактальной меры (1 (x) на единичном отрезке. С другой стороны, IFS отображений Tix = (1/2)x,T2x = (1/2)x + (1/2) с вероятностями pi,p2 позволяют генерировать инвариантную мультифрактальную меру ( (x) на аттракторе [0,1]. При этом (T1,p1), действуя на слово - точку x, генерируют символ (0) в префиксе, а (T2,p2) - символ (1). Последовательные итерации можно сделать марковскими, выбирая текущее Ti,i = 1, 2 в соответствии с матрицей переходных вероятностей. Полагая ( (x) ~ (1 (x), можно найти неизвестные вероятности и, следовательно, предсказать очередной символ в слове. Рассмотрим теперь подробно элементы описанной схемы.
1. Система итеративных функций
Пусть X = (X, й) - полное8 метрическое пространство, например, X с Кn. Отображение Б : X — X удовлетворяет условию Липшица, если для всех х,у £ X существует такое число с > 0, что
й (Б (х),) Б (у) < сй (х, у). (1)
Постоянной Липшица ЬгрБ называют минимальное с, для которого выполняется условие (1). Отображение Б называется сжатием, если ЬгрБ < 1. Принцип сжимающих отображений утверждает, что в полном метрическом пространстве (X, й) сжатие Б : X — X имеет единственную неподвижную точку хо = Б (хо) [32]. Более того, для любой точки уо £ X последовательность точек ук = Бок = Б(Бо(к-1)), к = 0,1, 2... сходится к хо в метрике й.
Пусть Н (X) - пространство непустых компактных подмножеств из X. Определим расстояние от точки х £ X до А еН (X) (рис. 1) как
й (х, А) = гп{уеАй (х, у). (2)
8То есть предел любой сходящейся последовательности Коши в X принадлежит X.
Тогда для любых A,B gH (X) выражение
dH (A, B) = sup {d (x, B), d (y,A) \x G A,y G B } (3)
называется метрикой Хаусдорфа и пара (H (X), dH) является полным метрическим пространством [21, 22, 33].
Конечный набор {иц}, г = 1 ,N сжимающих отображений иц : X —X с коэффициентами сжатия {ci}, i = 1,N называют системой гиперболических IFS [21, 23, 33-35]. Пусть Wi (A) = = {wi(x) \ix G A }. Коллективное действие IFS на A G H (X) описывает оператор Хатчинсона [21, 34, 35]
w (A) =1 w (A) (4) Рис. 1. Определение метрики Хаусдорфа. Точка
^-А=1 x G A наиболее удалена от B, точка y G B
Интуитивно понятно, что объединение
наиболее удалена от А
сжимающих отображений должно быть сжатием. Действительно, используя неравенство треугольника, легко показать [21, 22,35], что оператор w определяет сжатие
в (Н (X) )
dn (w (A), w (B)) < sdu (A, B) ,s = max {a} ,i = 1,N. (5)
Фундаментальная теорема теории IFS [21, 35] обобщает принцип сжимающих отображений на (H (X) , du ) :
Пусть {wi} ,i = 1,N система итеративных функций на X с Rn. Тогда существует единственное непустое инвариантное множество A G (H (X), dH) такое, что
A = w (A) = Wi (A) (6)
и для любой точки B G H последовательность B, w (B), w°2 (B),won (B) сходится к A в метрике Хаусдорфа dH при n — о. Множество A называют аттрактором IFS.
Уравнение (6) выражает свойство самоподобия аттрактора: A является объединением9 собственных, уменьшенных, копий
A = wi (A) U W2 (A) U ... U wn (A) = wi (wi (A) U W2 (A) U ... U wn (A)) U
UW2 (wi (A) U w2 (A) U ... U wn (A)) U ... U wn (wi (A) U w2 (A) U ... U wn (A))...
Большую роль в теории IFS играет теорема о коллаже [21, 32]: Пусть B G H (X) и {wi} ,i = 1,N - IFS, с максимальной константой Липшица c = max {ci}, имеет аттрактор A = w (A). Тогда
du (A, B) < (1 - c)"1 dH (B, w (B)). (7)
Иными словами, чем меньше расстояние (его называют коллаж-расстояние) между произвольным начальным множеством B и его образом w (B), тем ближе B к аттрактору A.
9Такое объединение называют коллажем.
2. Оператор Маркова: действие Т¥Б на меру
Представим себе массу, распределенную на некотором топологическом пространстве X. Мерой называют неотрицательное число, измеряющее ее количество в каждом подмножестве, из некоторого выделенного класса В (X) с X элементарных борелевых подмножеств [22, 36]. Этот класс содержит подмножества, которые можно представить комбинацией не более чем счетного числа объединений и пересечений открытых и замкнутых подмножеств из X. Примером для X С К служит набор интервалов вида (х^х^), [xi,xj), (х^х^], [х^х^], дополненный пустым множеством 0. Борелевой мерой ц на В (X) называется функция ц : В (X) — [0, то] такая, что ц (0) = 0 и для счетного набора борелевых непересекающихся множеств {Аг} £ В (X) , Ai П Aj = 0, г = ] выполняется условие аддитивности
411 г=1 АЧ =£ г=1 Ц(А).
Мера называется единичной или вероятностной, если ц (X) = 1. Частным случаем борелевой меры является мера Лебега - длина интервала. Однако мерой может быть не только длина. Борелеву меру можно ввести, используя динамику отображения / : [0,1] — [0,1]. Орбитой / длины п + 1 является последовательность точек хо,х1,х2,...,хп, где хп = /оп (хо) = / /о(п-1) (хо)). Пусть {1к }%=1, 1к = [(к — 1)/Х,к/Х) - разбиение единичного интервала. В результате итераций /оп (хо), некоторое число точек пк = 1к П хо,х1, ...,хп попадает в 1к. Если каждой итерации поставить в соответствие единицу времени, то пк - время пребывания орбиты в 1к. Можно надеяться, что при большом п и типичных начальных значениях хо числа рк (п, хо) = пк/(п + 1) не зависят от п и выбора хо. В этом случае почти для всех орбит их можно считать вероятностной мерой Цк для интервала 1к
рк (п,хо) & Цк = Иш (#(1к П{хо,х1,..,хп})/п), V Цк = 1.
п^те
Заметим, что точка у = / (х) принадлежит некоторому объединению интервалов и = 1т и 1Р и... и 13 тогда и только тогда, когда х £ /-1 (и). Это означает, что число точек последовательности х1, ...,хп+1, которые попали в и, совпадает с числом тех точек из хо,х1,..., хп, которые принадлежали /-1 (и). Числа Цк не должны зависеть от сдвига последовательности, так что выражение ц (и) = ц (/-1(и)) констатирует инвариантность меры.
Рассмотрим теперь вопрос о действии отображения на меру. В случае абсолютно непрерывной меры V существует интегрируемая функция (плотность) р такая, что V (I) = Л р (х) йх и для любой непрерывной функции ф интеграл относительно меры V определяется выражением
(ф, V) = J ф (х) р (х) йх.
Пусть Ь (х) кусочно-дифференцируемое отображение; А с I - интервал, который не содержит точек ' (х)| = 0, и Ь-1 (А) состоит из конечного числа интервалов Зк, каждый из которых монотонно отображается в А. Тогда для любой функции у = Ф (у) можно записать уравнение
Ф (y) dy = / ф (F (x)) |F' (x)| dx, ■JJk
J A JJk
где y = F (x). Полагая ф (y) = p (x) \1/F' (x)\, получим
/ p (x) |F' (x)| 1 dy = p (x) dx = v (Jk).
Jk
11 йу =
Если V = ц, то есть является инвариантной мерой, то действие на ц отображения Ь (суммирование по всем к ) дает уравнение Фробениуса - Перрона, определяющее плотность меры р (у),
р (у-) = Е Р {Хк=р(х) IVЬ' (х)|.
Для дискретного случая, который нас интересует, интеграл10 / gv от произвольной функции д (х) по мере V понимается как сумма
gv = V g (xk)m (xk) xk
m (xk ), (8)
где m (xk ) - мера, определенная в точке xk, и ряд сходится в абсолютном смысле [36]. Пусть {yi} значения функции y = f (xk). Их может быть несколько для каждой точки xk. Тогда меру точек yi естественно определить как меру их прообразов xk выражением
m (yi) = Xk=f-i{yi) m (xk )> (9)
которое имеет смысл, если число точек f-1 (yi) не бесконечно велико.
Символически действие (функции на меру называют push forward map и записывают как f * v = v (f, где правая часть является мерой точек x, которые отобразились в точки y = f (x). Для IFS с вероятностями {wf,pi} , i = 1, N, pi + .... + pn = 1 естественным обобщением такого действия является марковский оператор [21, 24, 35]
MW = Ei=1 PiV ◦ w-1, (10)
где (о) - знак композиции. Для любой непрерывной функции f : X ^ R справедливо выражение
fx fd(M(v)) = £N=1 Pi fx fd(v о w-1) =
Х* " ' ' (11)
= £ = рг/ы,(х) / (х) й(Ц ° Ш-1) =Т, = рг/х / ° (х)йЦ.
Заметим, что ц°ш~1 - мера тех точек х £ X, которые генерируют образы х — (х). Согласно (9), это можно записать как 8 (х — Шг (х)), что и приводит к замене / (х) — /ошг (х) в (11). Мера ц называется «самоподобной», если ее можно представить как линейную взвешенную комбинацию ее самой: ц = р1Ц°ш-1 + ...+рNц°ш—1.
10
Традиционная запись / д^ не имеет особого смысла для нашего случая.
3. Метрика Монжа - Канторовича в пространстве борелевых мер
Пусть М = М (В(X)) - пространство борелевых мер, и д : X ^ К функции с ЬХрд < 1. Метрика Монжа - Канторовича в М определяет расстояние между двумя мерами ^ и V как
dM (и, v) = sup
gdи - gdv
\\g (x) - g (y)\\<\x - y\, Vx,y e X}. (12)
Пара (M, dM) - полное метрическое пространство. Для любых и, v e M марковский оператор M является сжатием в (M, dM) [21,35]
dM (M(vi), M(v2)) < cdM (vi, V2).
Фундаментальная теорема утверждает [21, 35]:
Для IFS с вероятностями {wf, Pi},Yli=i Vi = 1 существует единственная бо-релева вероятностная мера и такая, что для любого подмножества A e B (X)
И (A) = M (A) = i V# о w-1 (A), i=i
и для любой меры vo e M последовательность vo, Mvo, Mo2vo,..., Monvo сходится к и в метрике dM, когда n — о. Носителем инвариантной меры является аттрактор A. Кроме того, для всех непрерывных функций g : X — R
J g (x) dи (x) = ^2N i Vi j g о Wi (x) dи (x). (13)
Численный вариант метрики Монжа - Канторовича известен как «экскаваторное расстояние» (earth-moving distance, EMD)n [37, 38]. Идея его вычисления заключается в следующем. Аппроксимируем две меры и и v зеркальными гистограммами с общим основанием. Пусть верхняя из них представляет собой кучи грунта и, насыпанного в каждый бин. Вторая, «опрокинутая» гистограмма, моделирует емкость ям v. Если площади гистограмм совпадают, то v способна вместить весь грунт и. Задача заключается в перемещении грунта из верхней гистограммы в ямы нижней с минимальными затратами. Последние зависят от выбранного «плана» перемещения и его «стоимости», которая определяется расстоянием dij (ground distance) между г-м бином, где находится грунт, и j-м бином ямы. Таким образом, вычисление EMD сводится к стандартной транспортной задаче линейного программирования (см., например, [39]) о минимизации транспортных расходов
Е™ i Е"= i fij • dij — min' Т!]= i xij = ai,
Е™ 1 !и = Ь], !ц > о, х = 1,т, 2 = 1 п,
где dij - стоимость перевозки единицы продукции из пункта X в пункт 2, - план перевозок X ^ 2, Ь] - потребности в продукте в пункте 2, ^ - запасы в пункте X. Для модели закрытого типа ^1 Ь] = ^™ 1 Основной трудностью решения задачи является большое число переменных и ограничений. Поэтому используются специальные алгоритмы [39, 40], которые с успехом адаптированы для вычисления расстояния между изображениями [41, 42].
11
Earth-moving machine - экскаватор.
4. Обратная задача IFS
Выше были описаны три ипостаси теоремы о сжимающем отображении. В евклидовом метрическом пространстве (X, d) сжатием является отображение Липшица (1), которое имеет единственную неподвижную точку. В пространстве компактов (H, du) сжатием в метрике Хаусдорфа du является оператор Хатчинсона (4) -объединение сжимающих отображений, а неподвижную точку называют аттрактором IFS. Наконец, в пространстве борелевых мер (M,dM) сжатие в метрике Мон-жа - Канторовича задает оператор Маркова (10) для IFS с вероятностями. Здесь неподвижной точкой является инвариантная единственная мера, носитель которой -аттрактор IFS.
Предположим, что мера V и аттрактор A известны. Обратная задача IFS состоит в нахождении параметров {wi} ,wi = cix + ai,i = 1,N и вероятностей {pi} или Pij. Рассмотрим вначале случай, когда выбор отображений в динамике случайных IFS производится независимо с помощью испытаний Бернулли, и опишем схему решения обратной задачи методом моментов, следуя работе [30]. Вопросы существования и сходимости решения подробно обсуждаются в статьях [43-45].
Для случайной переменной x G [0,1] определим набор статистических моментов
Vk = j xk dn, k = 0,1,2,..., (14)
где интеграл понимается в смысле формулы (8). Для вероятностной меры нулевой момент v0 = 1. В уравнении (13) выберем g (x) = xk. Тогда
(g ◦ wi) (x) = (Cix + ai)k = Vk J M ck-jxk-jak. (15)
Следовательно, уравнение (13) для инвариантной меры примет вид
^ ( к \ к-з к
Vk = j g(x) dv(x) = ^ N=1 k^)ck j ak Vk-j ■ (16)
Последнее выражение приводит к рекурсивной формуле для моментов
1 - E^ PiC0 Vk = E=o ( k ) Vk-j (EN, PiCk-jaj) ■ (")
Для марковской схемы приведенные формулы усложняются [53]. Прежде всего статистические моменты в (17) выражаются как сумма частичных моментов по числу IFS: Vn = Vn^ + ■■■ + Vn^ \ где Vj - решения системы уравнений, аналогичной (17),
Е (PjiCn - àij )VH) = - х0 ( k ) ck an^pj , (18)
i = 1,....,N,n > 1.
Здесь pji - марковские переходные вероятности и ôj - символ Кронекера. Нулевой момент для i-го IFS определяется выражением v0^ = mi, где mi - решения уравнений
, Pjimj = mi, i = 1,2,...,N, v0 = m, + .. + mN = 1. (19)
Теоретические моменты в (17)-(19) содержат искомые параметры: коэффициенты IFS и вероятности. Функционал обратной задачи основан на сравнении этих моментов с эмпирическими моментами щ, вычисленными по эмпирической мере (х, для к = 1, 2,.., M. В случае линейной оптимизационной задачи для функционала можно использовать Li-метрику (^|(ц — jij| ^ min), L2-метрику, коллаж-расстояние (7) или EMD. Задачу можно свести к квадратичной, используя «расстояние в шаре» (^M=1 (( k — (к)21к ^ min) [31]. Для численного решения используют стандартные оптимизационные пакеты или генетический алгоритм [46]. Критерием оптимального решения служит либо разность моментов (эмпирических и модельных), либо разности гистограмм эмпирической и модельной меры.
Для большинства практических задач проблема существенно упрощается. Так, например, если речь идет о предсказании пороговых значений некоторой величины, естественно использовать бинарный алфавит {0,1}. Следовательно, в уравнениях (17)-(19) достаточно ограничиться N = 2. Эмпирическая мера, полученная представлением слов в двоичном коде, имеет в качестве носителя единичный интервал. Поэтому можно выбрать IFS с фиксированными коэффициентами в виде w1(x) = (1/2)x, w2(x) = (1/2)x + 1/2 и аттрактором [0,1]. Обратная задача сводится тогда или к нахождению вероятностей pi, Р2, или в марковском варианте вероятностей pii и Р22 переходов символов 0 ^ 0 и 1 ^ 1, соответственно. Сопряженная пара вероятностей вычисляется из соотношений Р12 = 1 — Р11 и Р21 = 1 — Р22.
5. Практические приложения
Рассмотрим пример практического применения марковской IFS-схемы для прогноза магнитных бурь по вариациям индекса геомагнитной активности. Для тестирования модели предсказывались известные значения, которые не входили в обучающую выборку, то есть, фактически, делался эпигноз12.
Магнитные бури - это возмущения магнитного поля Земли, которые отслеживаются различными локальными и планетарными геомагнитными индексами [47]. Причинами бурь являются геоэффективные межпланетные возмущения солнечного происхождения: солнечные вспышки, выбросы корональных масс и межпланетные направленные взрывы, связанные с инверсиями межпланетного магнитного поля. Главная фаза геомагнитной бури представляет собой уменьшение Н-компоненты поля от 50 до 400 нанотеслов (нТ), которое продолжается от нескольких часов до суток и более. Депрессия вызвана пересоединением силовых линий межпланетного и геомагнитного полей и флуктуациями размеров магнитосферы под действием потоков солнечного ветра. Эти процессы приводят к ускорению имеющейся в ней плазмы до энергий порядка нескольких тысяч электрон-вольт, формируя таким образом кольцевой ток на расстоянии от 3 до 5 радиусов Земли. В северном полушарии в окрестности Земли он течет в западном направлении и порождает магнитное поле, направленное противоположно геомагнитному, и, следовательно, ослабляет его [48, 49].
В настоящей работе использован среднесуточный Д^-индекс13, который является мерой интенсивности кольцевого тока и вычисляется как усредненная величина
12Эпигнозом называют предсказание известных в прошлом значений, палегнозом (ро81^с1юп) -предсказание неизвестной истории ряда.
13Д8(;-индекс доступен на сайте http://swdcdb.kugi.kyoto-u.ac.jp/dstdir.
50 0
2 -50 а> ч
I -100
I
от
Р -150 -200 -250
Рис. 2. График среднесуточных значений Д^-индекса (1957-2002 годы). По вертикали, напряженность поля в нТ
0.03 т-
со
§■ 0.02 -
и:
со ^
о
(D Т
Q. 0.01 -s 1=
СО
0.0 -
0.0 0.06 0.12 0.18 0.24 0.31 0.37 0.43 0.49 0.55 0.61 0.67 0.73 0.79 0.85 0.92 0.98
Рис. 3. Гистограмма распределения бинарных слов K =12 для ряда Dst-индекса (1957-2002 годы) с порогом - 30 нТ
возмущений (рис. 2), отсчитываемых от спокойного уровня, по данным четырех магнитных обсерваторий, расположенных приблизительно вдоль магнитного экватора. В магнито-спокойные дни величина Dst лежит в пределах ± 20 нТ. Обычно, магнитной бурей считают значение Dst > —30 нТ; для сильных бурь Dst > —50 нТ и очень сильных14 Dst > —100 нТ.
Существуют различные подходы к предсказанию магнитных бурь [50-52]. В недавней работе [53] приведены результаты марковского прогноза магнитных бурь за период 1981-2002 годов по Dst-индексу на основе IFS и метода моментов. Ниже приводятся оригинальные результаты, полученные авторами для тех же данных двумя методами с целью проверки воспроизводимости результатов [53] и сравнения решений, полученных методами моментов и нашим методом, основанным на теореме о коллаже. Как отмечено выше, магнитная буря - сложный динамический процесс, однако ниже, как и в работах [18, 52], под «бурей» упрощенно понимается событие Dst > —30 нТ. На рис. 3 приведена оценка эмпирической меры, полученная по 16801 значениям временного ряда среднесуточных значений Dst-индекса для бинарного алфавита (порог - 30 нТ) и длины слова K=12, как и в [53]. Мульти-
14Во время сильной бури 11.02.1958 индекс Dst достиг величины - 409 нТ.
hluLilL-LuL
»* -' ■" Ш-u 11
0.8 . 0.6 -0.4 -0.2 ■ 0.0 -
0.0 0.2 0.4 0.6 0.8 1.0 а
Рис. 4. Мультифрактальные спектры больших отклонений для эмпирической меры -индекса: 1 - 1981-1991 годы, 2 - 1992-2002 годы
0.5 -0.4 -0.3 : 0.2 : 0.1 : 0.0 ;
0.0 0.2 0.4 0.6 0.8 1.0
Рис. 5. Кумулятивные гистограммы для эмпирической меры (кружки) и /^-модели (сплошная линия)
фрактальные спектры больших отклонений [54] для этой меры, полученные в пакете FracLab для двух различных интервалов времени (1981-1991 годы и 1992-2002 годы) практически совпадают (рис. 4). Таким образом, меру можно считать эмпирически стационарной.
Для построения модели использовался Д^-индекс с 1981 по 1996 год (всего 5844 значения). Обратная задача для IFS методом моментов решалась в среде MatLab 7.0 с использованием первых 15 моментов15 с помощью пакета Pattern search tool. Этот пакет содержит 5 методов поиска минимума оптимизационной задачи, включая генетический алгоритм. Все они дали оценки переходных вероятностей, совпадающие с точностью до 5 знаков после запятой. Наиболее устойчивый из полученных результатов16: рц = 0.91079904, р22 = 0.72421173. Для моделирования меры использовались 2 ■ 105 рекуррентных точек с указанными вероятностями.
В качестве второго метода решения обратной задачи использовалась теорема о коллаже в следующей модифицированной форме.
Оценка эмпирической меры |1 рассматривалась как неподвижная точка марковского оператора. Переходные вероятности находились минимизацией функционала d ((1, Mora (!) ^ min, где в качестве расстояния между двумя гистограммами H = {h%} и K = {к%} использовалась L\-метрика17 d (H, K) = Y1 % \h% — к%\. Полученные значения вероятностей рц = 0.91611578125, р22 = 0.69755109375 использовались для моделирования 2 ■ 105 рекуррентных точек модельной меры. Сравнение модели с эмпирической мерой проведено по графикам (рис. 5) соответствующих кумулятивных гистограмм18 F = £ % |h% — h|, построенных по аналогии с моделью случайного блуждания. Для тестирования модели использовался эпигноз 1095 значений нулей и единиц на временном интервале 1997-1999 годов. С практической точки зрения, ошибки в предсказании бури и ее отсутствия не равнозначны. Поэтому, для оценки качества прогноза использовались следующие три коэффициента, предложенные в работе [53]:
Т\ = п\/и, где Н\ - количество правильно предсказанных суффиксов19, п - общее количество прогнозов;
r2 = п2/п3, где п2 - количество правильно предсказанных суффиксов, если
15Именно это число моментов было использовано в [53].
16Напомним, что p12 = 1 — Pii, p21 = 1 — p22.
17Эксперименты показали, что EMD не дает здесь существенного преимущества.
18Здесь h - среднее по гистограмме.
"^-суффиксом называют предсказанные k = 1, 2,... символов в слове длиной K — k.
соответствующие реальные значения содержали единицу, п33 - количество прогнозов слов, реально содержащих единицу;
г3 = п4/п3, где п4 - количество «частично предсказанных» суффиксов, если соответствующие реальные значения содержали единицу. Частичное предсказание успешно, если независимо от длины суффикса была предсказана хотя бы одна единица. Таким образом, при вычислении и гз исключались слова, состоящие только из нулей. Полученные результаты для интервалов в 1-3 дня сведены в таблицу, где в качестве сравнения приведены значения коэффициентов для эпигноза по эмпирической мере. В последнем блоке жирным шрифтом приведены результаты, полученные в [53].
Таблица
Сравнение качества эпигноза на 1, 2 и 3 дня для Dst-индекса
Эмпирическая мера Коллаж-метод Метод моментов
% 1 день 2 дня 3 дня 1день 2 дня 3 дня 1 день 2 дня 3 дня
ri 86.43 76.62 69.29 86.98 77.26 69.66 86.52/75.83 77.26/58.12 69.66/43.73
Г2 58.00 27.75 14.61 70.00 35.08 19.63 68.67/67.24 34.03/34.89 19.18/19.47
Г3 58.00 50.79 46.12 70.00 59.69 54.79 68.67/67.24 59.16/55.78 54.79/48.44
Приведенные оценки показывают, что результаты вероятностного предсказания воспроизводимы и не слишком зависят от метода решения обратной задачи. Коллаж-метод дает лучшие результаты для краткосрочного предсказания.
Заключение
Марковское предсказание, основанное на решении обратной задачи IFS с вероятностями, эффективно в случаях, когда в отсчетах временного ряда можно обоснованно выделить несколько диапазонов величин, согласованных с некоторым «грубым» разбиением пространства состояний динамической системы. Процедура марковского предсказания содержит следующие этапы.
• Выбор объема |2| алфавита, порога (уровней) и длины слова K. Порог, или уровни, выбираются обычно из физических соображений; их число определяет |2|. Длина слова K должна обеспечивать статистически значимое заполнение всех бин гистограммы.
• Построение эмпирической меры. Проверка ее стационарности и наличия мультифрактальных свойств. Стационарность меры проверяется совпадением гистограмм, построенных для двух фрагментов ряда. Мультифрактальные свойства меры обнаруживаются вычислением лежандровского спектра, либо спектра больших отклонений.
• Решение обратной задачи для IFS с вероятностями на основе эмпирической меры. Число отображений N, входящих в IFS, выбирается равным объему алфавита N = |2|. При фиксированных коэффициентах IFS в методе моментов свободным параметром является номер M максимального момента. Увеличение M уменьшает ошибку функционала. Однако чрезвычайно большие M могут приводить к неустойчивости численного алгоритма, поскольку xM, x < 1 умножаются в (17), (18) на большие значения биноминальных коэффициентов. Результатом решения является матрица переходных вероятностей.
• Генерация модельной меры с помощью рекуррентной последовательности IFS с полученными значениями pгj с объемом, превышающим объем эмпирической меры. Проверка модели реализуется сравнением гистограмм, эмпирической и модельной.
• Модельная мера, на которой основано предсказание, должна хорошо представлять все возможные слова. В случае бинарного алфавита, успешное предсказание получается, когда гистограмма имеет форму «корыта», приводящего к сравнимым вероятностям pll и p22. Для прогноза на один символ20 рассматривается фрагмент слова длины K — 1 : sl,s2, ...sk-l. Альтернатива (0 или 1) для sk выбирается на основе модельной меры как max{P(si, s2, ...sk-i, 1), P(si, s2, ...sk-l, 0)}.
Обобщение на прогноз двух и более символов очевидно. В качестве возможного варианта можно делать упрощенный прогноз, опираясь только на эмпирическую меру. Однако его качество всегда хуже, поскольку в реальной выборке могут быть плохо представлены или вообще отсутствуют некоторые из возможных слов.
Приведенный практический пример показывает, что коллективная линейная динамика гиперболических отображений с вероятностями может быть полезна для обоснованного вероятностного предсказания сложных природных процессов. Возможные использования обобщенных фрактальных преобразований [55] и IFS с коэффициентами сг > 1/2 еще ждут своих исследователей.
Библиографический список
1. Noakes L. The Takens embedding theorem II Inter. J. Bifurcation and Chaos. 1991. Vol. 1. P. 8б7.
2. Sauer T., YorkeJ.A., CasdagliM. Embedology Il J. Statist. Phys.1991. Vol. б5. P. 579.
3. Афраймович B.C., Рейман A.M. Размерности и энтропии в многомерных системах ИНелинейные волны. Динамика и эволюция. М.: Наука, 1989. C. 238.
4. Макаренко Н.Г. Реконструкция динамических систем по хаотическим временным рядам II Нелинейные волны'2004. Нижний Новгород, 2004. С. 398.
5. Stark /.Delay reconstruction: dynamics versus statistics II Nonlinear dynamics and statistics IA.I. Mees ed. Birkhauser, 2001. P. 81.
6. Макаренко Н.Г. Эмбедология и нейропрогноз II Лекции по нейроинформатике. Ч. 1. Нейроинформатика-2003. Москва, 2003. C. 8б.
7. Poggio T., Girosi F. A theory of networks for approximation and learning Il MIT AI Lab. Techn. Rep. 1989. Memo №1140. https^Ihpdsl.mitedulbitstreamlmU^l 1I2IAIM-1 l40.pdf
8. Малинецкий Г.Г., Потапов А.Б. Современные проблемы нелинейной динамики. М.: УРСС, 2002. 358 с.
9. Farmer J.D., Sidorovich /./.Predicting chaotic time series II Phys. Rev. Lett. 1987. Vol. 59. P. 845.
l0. Kantz H., Schreiber Th. Nonlinear time series analysis. Cambridge Univ.Press, 2004. 3б9 p.
20В случае бинарного алфавита.
11. McSharry P.E. Innovations in consistent nonlinear deterministic prediction. D.Phil. Thesis. University of Oxford, 1999.
12. Nakamura T., Kilminster D., Judd K., Mees A. A comparative study of model selection methods for nonlinear time series // Int. J. of BifUr. and Chaos. 2004. Vol. 14. P. 1129.
13. Mukhin D.N., Feigin AM.,Loscutov E.M., Molkov Y.I. Modified Bayesian approach for the reconstruction of dynamical systems from time series // Phys.Rev.E. 2006. Vol. 73(3 Pt 2):036211.
14. Kantz H., Ragwitz M. Phase space reconstruction and nonlinear predictions for stationary and nonstationary Markovian processes // Intern. Journal of Bifurcation and Chaos. 2004. Vol. 14, № 6. P. 1935.
15. Froyland G. Extracting dynamical behaviour via Markov models //Nonlinear dynamics and statistics / A.I. Mees ed. Birkhauser, 2001. P. 283.
16. Froyland G. Markov modelling for random dynamical systems. 1998. http://www.maths.unsw.edu.au/froyland
17. Daw C.S., Finney C.E.A., Tracy E.R. A review of symbolic analysis of experimental data // Rev. of Scientific Instruments. 2003. Vol. 74. P. 916.
18. Wanliss J.A., Ahn V.V., Yu Z.G., Watson S. Multifractal modeling of magnetic storms via symbolic dynamics analysis // J. Geopys. Res. 2005. Vol.110. P. AO814.
19. Anh Vo, Lau Ka-Sing, Yu Zu-Gao. Multifractal characterization of complete genomes // J.Phys. A: Math.Gen. 2001. Vol. 34. P. 7127.
20. Tino P. Multifractal properties of Hao's geometric representations of sequences // Physica A. 2002. Vol. 304(3-4). P. 480.
21. Barnsley M. Fractals everywhere. N.Y.: Academic Press, 1988. 531p.
22. Falconer K. Fractal geometry. Mathematical Foundations and Applications. Wiley, 2003. 337 p.
23. Barnsley M.F., Demko S. Iterated function systems and the global construction of fractals // Proc. Roy. Soc. London A. 1985.Vol. 399. P. 243.
24. Hutchinson J.E. Fractals and self-similarity // Indiana Univ. Math. 1981. Vol. 30. P. 713.
25. Diaconis P. Iterated random function // SIAM Review. 1999. Vol. 41. P. 45. http://www.stat.berkeley.edu/census/511.pdf
26. Vrscay E.R. From fractal image compression to fractal-based methods in mathematics // Fractals in Multimedia /ed. by M.F. Barnsley, D. Saupe and E.R. Vrscay. New York: Springer-Verlag, 2002.
27. Jeffrey H.J.Chaos game representation of gene structure// Nucleic Acids Research. 1990. Vol. 18. P. 2163.
28. Tino P. Spatial representation of symbolic sequences through iterative function systems // IEEE Transactions on Systems, Man, and Cybernetics: Part A: Systems and Humans, 1999. Vol. 29(4). P. 386.
29. Tino P., Dorffner G. Predicting the future of discrete sequences from fractal representations of the past //Machine Learning. 2001. Vol. 45(2). P. 187. http://www.cs.bham.ac.uk/pxt/my.publ.html.
30. Barnsley M.F., Ervin V., Hardin D., Lancaster J. Solution of an inverse problem for fractals and other sets //Proc. Nat. Acad. Sci. USA. 1986. Vol. 83. P. 1975.
31. Iacus St. M., Torre D.L. Approximating distribution functions by iterated function systems// Departemental Working Papers 2002-03, Department of Economics University of Milan Italy. http://ideas.repec.org/e/pla155.html.
32. Hart J.C. Computer display of linear fractal surfaces // Doctor Thesis. University of Illinois at Chicago, 1991. http://graphics.cs.uiuc.edu/jch/papers/diss.pdf.
33. Макаренко Н.Г. Фракталы, мультифрактальные меры и аттаракторы // Нелинейные волны'2002. Нижний Новгород, 2003. С. 381.
34. Макаренко Н.Г. Фракталы, аттракторы, нейронные сети и все такое // Лекции по нейроинформатике. Ч.2. // Нейроинформатика-2002. Москва, 2002. С. 121.
35. Falconer K. Techniques in fractal geometry. Wiley & Sons, 1997. 256 p.
36. Hutchinson J.E. Measure Theory. 1995. http://wwwmaths.anu.edu.au/john/lecture_notes.html.
37. Rubner Y., Tomasi C., Guibas L.J. The Earth mover's distance as a metric for image retrieval // Technical Rep. STAN-CS-TN-98-86.
38. Kaijser T. Computing the Kantorovich distance for images // J. Mathematical Imaging and Vision. 1998. Vol. 9. P. 173.
39. Лемешко Б.Ю. Методы оптимизации. Конспект лекций. http://www.ami.nstu.ru/ headrd/
40. Stark J. A neural network to compute the Hutchinson metric in fractal image processing // IEEE Trans. Neural Networks. 1991. Vol 2. P. 156.
41. Wadstromer ^.Coding of fractal binary images with contractive set mappings composed of affine transformations // PhD Theses. Linkopings univer. 2001.
42. Ling H., Okada K. EMD-L1: An efficient and robust algorithm for comparing histogram-based descriptors //European Conference on Computer Vision. 2006. http://www.cs.umd.edu/ hbling/main.htm.
43. Forte B., Vrscay E. R. Solving the inverse problem for function/image approximations using iterated function systems. I.Theoretical basis; II. Algorithm and computations // Fractals. 1994. Vol. 2,3. P. 325; P. 346.
44. Handy C.R., Mantica G. Inverse problems in fractal construction: moment method solution // Phys. D. 1990. Vol. 43. P. 17.
45. Abendat S., Demko S., Turchetti G. Local moments and inverse problem for fractal measures // Inverse Problems. 1992. Vol. 8. P. 739.
46. Lutton E., Levy-Vehel J., Cretin G., Glevarec Ph., Roll C. Mixed IFS: Resolution of the inverse problem using genetic programming // Complex Systems. 1995. Vol. 9. P. 375.
47. Заболотная Н.А. Индексы геомагнитной активности. М.: Гидрометиздат, 1977. 39 с.
48. Яновский Б.М. Земной магнетизм. Ленинград: ЛГУ, 1978. 592 с.
49. Пудовкин М.И., Распопов О.М., Клейменова Н.Т. Возмущения электромагнитного поля Земли. Ленинград: ЛГУ, 1976. 247 с.
50. Watanabe Sh., Sagawa E., Ohtaka K., Shimazu H. Prediction of the Dst index from
solar wind parameters by a neural network method // Earth Planets Space. 2002. Vol. 54. P. 1263.
51. Stepanova M., Antonova E., Troshichev O. Prediction of Dst variations from Polar Cap indices using time-delay neural network // J.Atmosph. and Solar-Terrestrial Phys. 2005. Vol. 67. P. 1658.
52. Strivastava N. A logistic regression model for predicting the occurrence of intense geomagnetic storms //Ann.Geophys. 2005. Vol.23. P.2969.
53. Ahn V.V., Yu Z.G., Wanliss J.A., Watson S.M. Prediction of magnetic storm events using the Dst index // Nonlinear Processes in Geophysics. 2005. Vol. 12. P. 799.
54. Levy Vehel J. Numerical computation of the large deviation multifractal spectrum // URL: http://www-rocq.inria.fr/fractales.
55. Forte B, Vrscay E.R. Theory of generalized fractal transforms. Fractal Image Encoding and Analysis / Edited by Y. Fisher. Heidelberg: Springer Verlag, 1998. http://links.uwaterloo.ca/person.ed.html
Институт математики, Алма-Ата, Поступила в редакцию 19.06.06
Казахстан После доработки 15.07.06
Санкт-Петербургский государственный
университет
ITERATED FUNCTION SYSTEM AND MARCOVIAN PREDICTION OF TIME SERIES
N.G. Makarenko, L.M. Karimova, S.A. Mukhamejanova, I.S. Knyazeva
This paper demonstrates a tool for prediction time series on a base of iterated function system of the theory of fractals. Iterations result in an attractor or fractal in a space of compacts. The attractor is a support of invariant probabilistic measure or multifractal in a space of Borel measures. An inverse problem consists of finding iterated function system and its probabilities by means of empirical measure. The estimates might be obtained from time series by symbolic dynamics methods. In addition to necessary mathematical material two practical results of predictions of threshold values for financial time series and geomagnetic storms are represented.
Макаренко Николай Григорьевич - родился в Троицке (1945). Окончил Уральский государственный университет им. Горького (Свердловск) по специальности «астрономия». Ведущий научный сотрудник, д.ф.-м.н. ГАО РАН (Пулково, Санкт-Петербург). Главный научный сотрудник, д.т.н. Лаборатории компьютерного моделирования (Институт математики, Алма-Ата, Казахстан). Область научных интересов: фрактальная геометрия, вычислительная топология, детерминированный хаос, нейронные сети, физика Солнца. Имеет более 70 научных публикаций. Е-таП:^-такаг@таП.ги
Каримова Лаиля Митхатовна - родилась в Алма-Ате (1948). Окончила Ленинградский государственный университет по специальности «радиофизика». Старший научный сотрудник, к. ф.-м.-н., Лаборатории Компьютерного Моделирования (Институт Математики, Алма-Ата, Казахстан) Область научных интересов фрактальная геометрия, нейронные сети, математическая морфология. Имеет более 40 научных публикаций.
Мухамеджанова Светлана Адиковна - 1979 года рождения, окончила механико-математический факультет Казахского Государственного Национального Университета, младший научный сотрудник Лаборатории Компьютерного Моделирования (Институт Математики, Алма-Ата, Казахстан) Область научных интересов - фрактальная геометрия, нейронные сети, прогноз временных рядов.
Князева Ирина Сергеевна - 1983 года рождения, студентка 2-го курса магистратуры физического факультета Санкт-Петербургского государственного университета. Область научных интересов - фрактальная геометрия, нейронные сети, прогноз временных рядов.
