CC BY
37
ВЫЯВЛЕНИЕ НЕЛИНЕЙНЫХ СВЯЗЕЙ МЕЖДУ СТОХАСТИЧЕСКИМИ ОСЦИЛЛЯТОРАМИ ПО ВРЕМЕННЫМ РЯДАМ
Д.А. Смирнов
Рассматривается задача выявления и количественной оценки нелинейных направленных связей между стохастическими осцилляторами. Предложены характеристики связи и методика их оценки по временным рядам. Аналитически получено выражение для уровня статистической значимости вывода о наличии связи, что позволяет делать надежные заключения при анализе относительно коротких сигналов. Эффективность подхода демонстрируется в численных экспериментах на примерах осцилляторов с различными свойствами и различными видами функций связи.
Ключевые слова: Оценка связей, связанные осцилляторы, нелинейный анализ временных рядов, фазовая динамика, статистические выводы.
Введение
При анализе сложных систем различной природы очень важен вопрос о наличии и характере взаимодействий (связей) между их подсистемами. В частности, структура и интенсивность связей в ансамбле колебательных систем определяет возможность их синхронизации, что является предметом многочисленных исследований (см., например, [1-5]). В последние годы активно исследуются возможности выявления направленных связей по временным рядам [6-16], что актуально при исследовании колебаний в сердечно-сосудистой системе [8,17-20] и почках [21,22], процессов переноса информации в нейронных сетях [23-30] и пр.
Один из наиболее чувствительных методов оценки связи между двумя осцилляторами был предложен в [6]. Он основан на моделировании фазовой динамики и применим в случае достаточно длинных временных рядов (несколько сотен характерных периодов) или малых шумов. В работе [7] метод обобщен на случай более коротких временных рядов и существенных шумов за счет специальных поправок и эмпирически найденного порогового значения характеристики связи, соответствующего 95%-й доверительной вероятности вывода о наличии связи. Метод нашел применения в нейрофизиологии [31] и климатологии [32]. Он был далее развит в
работе [33], где предложены модифицированные характеристики связи и получено аналитическое выражение для доверительной вероятности, с которой можно сделать вывод о наличии связи. Однако оба упомянутых способа улучшения метода ориентированы на конкретную форму используемой эмпирической модели, а именно, на связи, описываемые членами невысокого порядка в уравнениях фазовой динамики. Они не пригодны для выявления нелинейных воздействий произвольной сложности.
Выявление нелинейности зачастую важно для того, чтобы понять физический механизм взаимодействия. Кроме того, нелинейность связей может существенно влиять на динамику [34]. Поэтому в данной работе развитый в [33] метод (п. 1.1) обобщается для оценки нелинейных связей любого порядка (п. 1.2). Применимость новых оценок связи в случае достаточно коротких временных рядов и осцилляторов с различными свойствами иллюстрируется в численном эксперименте (раздел 2). Выводы представлены в Заключении.
1. Количественные характеристики направленных связей и их оценки
Основная идея подхода - оценить, насколько сильно зависит будущая эволюция фазы одной системы от текущего значения фазы другой системы. Формализм представлен ниже.
1.1. Характеристики первого порядка и их оценки
1.1.1. Теоретические характеристики. Фазовая динамика слабо связанных автоколебательных систем (детерминированных, в отсутствие связи периодических) с хорошей степенью точности [35] описывается системой дифференциальных уравнений вида
= Ю1 + а1 8ш(ф2 — Ф1), = Ю2 + а1 8т(ф1 — Ф2), где юи - собственные частоты осцилляторов, аи - коэффициенты связи (к = 1, 2). Если осцилляторы слабо возмущаются шумом, то модель может быть обобщена [36] и записана в виде системы стохастических дифференциальных уравнений
(ф1 /(И = Ю1 + а1 8т(ф2 — Ф1) + (¿),
(2)
(ф2/(М = Ю2 + а1 8т(ф1 — Ф2) + ^2
где ^ (¿) - независимые источники белого шума с нулевым средним и автоковариационными функциями (^(£')) = о^— ¿'), параметр о^ характеризует интенсивность шума.
Пусть величины о^к и \аи\ достаточно малы, так что вклад соответствующих слагаемых в уравнении (2) в приращение фазы Дфи (¿) = фи (£ + т) — фи (¿) мал по сравнению с «линейным приращением» шд.т, где т - конечный интервал порядка характерного периода колебаний или больше. Тогда можно от уравнений (2) перейти к разностным уравнениям, что требуется при дальнейшей работе с дискретными временными рядами. Этот переход достигается путем интегрирования уравнений (2) на интервале длиной т, что дает
Дф1^) = W1 + а1 ОО8(ф2С0 — ф1С0) + в1 В1п(ф2(*) — ф1 СО) + £1(г),
Дф2^) = W2 + а2 СО8(ф1^) — ф2(^)) + в2 81п(фх(^) — ф2(¿)) + £2®,
где £k(t) - шумы с нулевым средним и дисперсиями о^ = о|кт. Их автокорреляционные функции (АКФ) у^(l) = (^k(t)efc(t + l))/o"lk линейно спадают от единицы до нуля на интервале временных лагов 0 < l < т и равны нулю при l > т.
Эти свойства следуют из того, что при малости шумов и связи имеет место [7,33]
t+т
£к(t) — j Ik(t')dt'.
t
Сила воздействия j-го осциллятора на k-й (j ^ к) определена в [33] из следующих соображений. Рассмотрим статистические свойства приращения фазы Дфk -левой части уравнений (3). Ее среднее значение (Дфk) = Wk — ®kт- Из-за влияния шума и другого осциллятора величина Дфk флуктуирует около среднего значения (модуляция периода колебаний). При слабой связи стационарное вероятностное распределение «свернутых» фаз (ф1 mod 2п, ф2 mod 2п) является практически равномерным в квадрате [0, 2п) х [0, 2п). Значит, функции-слагаемые в правой части (3) (это 1, cos^j — фk) и sin^j — фk)) взаимно ортогональны в этой области, то есть (cos^j — фk) sin^j — фk)) = 0 и т.д., где угловые скобки означают математическое ожидание (усреднение по указанной равномерной вероятностной мере). Возведение в квадрат обеих частей уравнения (3) и усреднение показывает, что дисперсия Дфk -это сумма двух слагаемых
оДфк = Cj^k + о^, (4)
где величина
Cj^k = 2 (afc + ßk) (5)
названа силой воздействия j ^ к.
Таким образом, согласно выражениям (4) и (5) интенсивность флуктуаций приращения фазы Дфk представлена в виде суммы двух независимых факторов - влияния шума и другого осциллятора. Это дает ясный физический смысл величине Cj^k.
1.1.2. Статистические оценки. На практике уравнения (3) неизвестны, а имеются временные ряды от двух систем - {x1(t1), ...,x1(tNx)} и {x2(t1), ...,x2(tNx)}, где Xk - наблюдаемые, ti = iДt, Дt - интервал выборки, Nx - длина ряда. Чтобы оценить силу связи Cj^k, рассчитывают фазы колебаний {ф^^),..., ф^^)} и {ф2(t1), ф2(tN)} с помощью известных методов [37,38] и оценивают коэффициенты модели (3) путем минимизации среднего квадрата ошибки
S (wk, ak, ßk) =
1 N-т/дt 2
= 1Тт-7X7 ^ (Дфk(ti) — Wk — ak cos^j (ti) — фk(ti)) —ßk sin(фj(ti) —фk(ti))) .
N —т/Д i=1
Таким образом, оценки коэффициентов Wk, ak, ßk есть arg min S(wk, ak, ßk). В ка-
Wk >ak>ßk
честве оценки величины Cj^k можно использовать Cj^k = (1/2) (ak + ß^.
Если временной ряд очень длинный, то Wk, ak, ßk практически совпадают с истинными значениями Wk, ak, ßk, а значит, и Cj^k = Cj^k. Однако на практике чаще встречаются ряды весьма умеренной длины. При этом Wk, ak, ßk уже не равны в точности истинным значениям коэффициентов, так что необходимо указать вероятную величину погрешности оценок. Уровень статистической значимости (вероятность случайной ошибки) вывода о наличии воздействия j ^ k, то есть вы-
вода о том, что с-^и > 0, рассчитывается следующим образом. При гипотезе отсутствия связи (с-^и = 0) случайные величины аи и |3д. являются независимыми и одинаково распределенными по нормальному закону с нулевым средним и некоторой дисперсией о0_к. Следовательно, величина Х^и = (а| + Ри )/о0.к распределена при этом по закону «хи-квадрат» с двумя степенями свободы. Ее функция распределения Ф2(х) табулирована (см., например, [39]). Обозначим Х2,1—р такое число, что Ф2(12>1—р) = 1 — р, то есть (1 — р)-квантиль распределения. Если оказывается, что Х^и > I2 1—р, то можно опровергнуть гипотезу об отсутствии связи на уровне значимости р. Чем меньше р, тем надежнее вывод. Обычно значение р = 0.05 счита-
О2
но вместо нее можно использовать оценку, полученную в [7],
ется соответствующим достаточно высокой надежности. Величина оа неизвестна,
ак
• 2 20вк Л l \ ( Wk + Wj) l \ ( l (С +
0а = ^ (l + 2E 1 - 7/At СОП т/Atj exp '
ак
N \ т/At) V т/At / l 2т/At
где оценка дисперсии шума о2 = min S(wk, , ßk). Таким образом, уровень зна-
wk,ak ,ßk
чимости, на котором можно сделать вывод о наличии воздействия j ^ к, оценивается как pj^k = 1 — Ф2(х,2^). Формула для оценки о|к получена при предположении о том, что АКФ шума £k линейно спадает до нуля на интервале [0, т]. Это единственный момент, где используется предположение о свойствах шума, но он важен, так как иначе нельзя оценить погрешности оценок связи.
1.2. Характеристики более высокого порядка и их оценки
1.2.1. Теоретические характеристики. Основное внимание в данной работе направлено на случай нелинейностей более высокого порядка в уравнениях фазовой динамики. А именно, функции связи в уравнениях (1) и (2) могут зависеть не только от разности фаз из-за сильной нелинейности исследуемых систем и их взаимных воздействий. В последнем случае необходимо рассматривать более общую модель фазовой динамики [37]
(ф1М = Ю1 + С1(фь ф2)+ Ы^,
(6)
(ф2М = «2 + С2(ф2, ф1) +
где функции Си 2п-периодичны по обоим аргументам. Они определяют взаимодействие осцилляторов и собственную нелинейность их динамики. Как и в предыдущем случае, можно перейти к следующим разностным уравнениям при условии малости шумов ^и(1) и величин \Си\:
Дф1^) = ф2(^), а1) + £1(1),
(7)
Дф2(^) = ^2(ф2(^), ф1(г), а2) + £2(г),
где Еи - тригонометрические многочлены вида
Р:(фи, ф-, аи) = Wk + ^ (аи,т,п СО8(тфи — пф-) + ви,т,п 81п(шфи — пф-)) ,
где к = 1,2; ak = (wk, {afc
)eQk) - вектоРы к°эффщиентов; Qk - диапазон суммирования, то есть набор пар значений m и n, определяющих, какие слагаемые присутствуют в многочлене. При близких частотах колебаний осцилляторов члены с m = n = 1 являются следствием линейной связи вида или к2^1(x2 — x1). Члены с n = 2 могут отражать воздействие, квадратичное по координате влияющего осциллятора, например, k2^\x2- Возможны и различные комбинации, так что связь в уравнениях фазовой динамики может описываться набором одночленов различных порядков с n = 0. Наибольшее влияние на динамику оказывают резонансные члены, то есть соответствующие m/n ~ Wjв уравнении для k-го осциллятора. Однако нерезонасные слагаемые тоже могут быть важны [34].
Ниже предлагается определить силу воздействия j ^ к по аналогии с подходом [33] и выводится аналитическая формула для уровня значимости в рассматриваемом нелинейном случае. Рассмотрим вновь статистические свойства приращения фазы Дфд. (в левой части уравнений (7)) при условиях слабой связи между осцилляторами, слабой нелинейности и малого шума. Среднее значение (Дф&) = wk ~ Стационарное вероятностное распределение «свернутых» фаз (ф1 mod 2п, ф2 mod 2п) является вновь равномерным в квадрате [0, 2п) х [0, 2п), а слагаемые в многочлене Fk при равномерной плотности распределения являются взаимно ортогональными функциями в этой области. Поэтому возведение в квадрат обеих частей (7) и усреднение дает
(O^k )2) = w2 + 1 £ (al,m,n + вк,т,п) + < • (9)
(m,ri)€Qk
Слагаемые в (9), соответствующие n = 0, описывают воздействие j-го осциллятора на k-ый. Обозначим их сумму Oj^k и назовем ее силой воздействия j ^ к
Cj^k = 77
Z (ak
,m,n k,m,n
(10)
(m,n)£Qk, n=0
Остальные слагаемые определяют нелинейность индивидуальной фазовой динамики. Обозначим их bk
bk =2 Е (ak,m,n + $k,m,n)^ (11)
(m ,n)€Qk, n=0
Поскольку (Дфk) = wk, дисперсия приращения фазы оДфк = ^(Дфд;)2^ — (Дфk)2 определяется следующими тремя слагаемыми:
оДфк = bk + Cj^k + о^ • (12)
Величину Cj^k можно по-разному нормировать. Во-первых, Cj^k/w2 показывает силу воздействия по сравнению с линейным приращением фазы, то есть интенсивность модуляции периода колебаний за счет воздействия j-го осциллятора. Эта величина всегда много меньше единицы при условии малости связи. Во-вторых, Cj^k/оДфк показывает долю флуктуаций величины Дфk, обусловленную воздействием j-го осциллятора. Эта величина меньше единицы, но может почти достигать ее, если фазовая
нелинейность и шум малы по сравнению с оцениваемым воздействием. В-третьих, cj^k/°efc показывает силу воздействия j-го осциллятора по сравнению с воздействием шума.
Введенная характеристика связи зависит от параметра т - временного масштаба. На различных временных масштабах относительные роли шума и воздействия j-го осциллятора могут быть различны. Так, влияние белого шума или шума, время автокорреляции которого много меньше характерного периода колебаний осцилляторов T, преобладает при т ^ T .С ростом т растет роль остальных слагаемых в (12), но при очень больших т вновь доминирует шум, так как зависимость приращения фазы от ее текущих значений пропадает. Таким образом, при промежуточных значениях т влияние одного осциллятора на другой проявляется наиболее сильно. Практика показывает, что целесообразно брать в качестве т величину основного периода колебаний T [6,7]. Если собственные периоды колебаний двух осцилляторов Ti и T2 сильно различаются, то можно рассматривать характеристики связи при двух значениях т. В разделе 2 для краткости приводятся результаты только для одного значения т: т = (Ti + T2)/2 в случае близких Ti и T2 или т = min{Ti,T2}, если Ti и T2 сильно различаются.
1.2.2. Статистические оценки. Чтобы реализовать предложенный формализм при анализе временного ряда, зададим достаточно большой набор одночленов в (8), чтобы он содержал все слагаемые, необходимые для описания исследуемого взаимодействия. Согласно методу наименьших квадратов минимизируем
N-т/At
^ (ak ) = N -т/м £ (ti) - Fk (^ (ti), Ф (ti), ak ))2. i=i
Тогда a2, = min S (ak) и äk = arg min S (ak) - оценки дисперсии шума и коэф-
k ak ak
фициентов многочлена Fk, соответственно. Оценки характеристик связи получим аналогично (5)
cj^k = 2 У ] i^ßk,m,n + ßk,m,n) • (13)
(m,ra)GQk, n=0
Как и в п. 1.1, величина Cj^k всегда неотрицательна. Даже при отсутствии связи между осцилляторами, то есть Cj^k = 0, оценка ßj^k почти наверное принимает значения, большие нуля. Таким образом, достоверный вывод о наличии воздействия j ^ k нельзя сделать из условия Cj^k > 0. Его можно сделать, только если Cj^k «существенно» больше нуля. Для получения количественного критерия, найдем закон распределения Cj^k при отсутствии связи и нелинейности, то есть для системы (6) с Gi = G2 = 0. В этом случае оценки всех коэффициентов äk,m,n, ß k,m,n статистически независимы друг от друга и распределены одинаково по нормальному закону с нулевым средним [7]. Их дисперсии обозначим a|k . Учитывая, что сумма квадратов M независимых нормально распределенных с нулевым средним и единичной дисперсией величин распределена по закону «хи-квадрат» с M степеня-
а2 +ß2
ми свободы, получим, что величина = Y1 k'mJi—^^ распределена по
m,n(n=0) ak,m,n
закону «хи-квадрат» с Mk степенями свободы, где Mk - количество слагаемых в
выражении (13). Значение и можно рассчитать по временному ряду, если вместо
дисперсии o|fc подставить оценку [7]
а2 fl+2 V U- —) cos (l{m™k+j \ exp ( ^^ ^]
— N 11+2 ¿i I1 т/AiJ т/Ai ) p I 2x/At
Таким образом, получен критерий для того, чтобы судить о статистической значимости вывода о наличии воздействия j — k. А именно, если величина Х2^, рассчитанная по временному ряду, превосходит (1 — p) -квантиль распределения «хи-квадрат» с Mk степенями свободы, то вывод можно сделать на уровне значимости p, то есть вероятность случайного ошибочного обнаружения связи не превосходит p. Функция распределения для закона «хи-квадрат» с Mk степенями свободы Фмк табулирована [39]. Зная ее, можно оценить уровень значимости, на котором X^k отлична от нуля, по формуле pj^k - 1 — Фмк (Xj^k)•
Сравнивать силы воздействия j — k при различных условиях непосредственно по величинам Cj^k (13) не удается, так как последние являются смещенными оценками истинных значений Oj^k (10), что можно показать аналогично [7], то есть (Cj^k) > cj^k, где угловые скобки означают математическое ожидание. Чтобы исключить это смещение (систематическую ошибку), достаточно вычесть из квадрата оценки каждого коэффициента в (13) ее дисперсию аналогично [7]. Полученную несмещенную оценку обозначим
Cj^k = 2 ^ (ak,m,n + bk,m,n — 2alk,m,n) • (14)
m,n (n=0)
Суммируя, отмечу, что предложенные оценки и формула для уровня значимости применимы при следующих условиях.
1) Хорошо определены фазы колебаний. Это имеет место в случае узкополосного спектра мощности сигналов, когда можно применять любой из известных методов расчета фазы, например, с использованием преобразования Гильберта [37]. Для более сложных сигналов введение фаз должно быть обеспечено с помощью специальных подходов. В примерах раздела 2 фазы хорошо определены и эта проблема подробно не обсуждается.
2) Слабая зависимость между фазами осцилляторов, то есть почти равномерное распределение «свернутой» разности фаз на отрезке (0, 2п). Это условие можно контролировать, рассчитывая по временному ряду обобщенные коэффициенты фазовой синхронизации pm,n — |(et(m^l(i)_n^(i))^>J, где m,n - положительные целые числа. Предложенные оценки связи применимы для данного ряда, если pm,n < pc для всех m, n из множества Qk. Если pm,n > pc для некоторой пары индексов m, n, то временной ряд не пригоден для оценки направленных связей. Если рассматриваются осцилляторы с близкими собственными частотами, то для них большое значение может принимать только pi,i. Эта величина обозначается далее просто p. Пороговое значение pc — 0-45 подобрано эмпирически на эталонных примерах (см. раздел 2).
3) Шум £k имеет АКФ, линейно спадающую от единицы до нуля на интервале временных лагов [0, т] и равную нулю при больших лагах. Для проверки следует
оценить АКФ остаточных ошибок - слагаемых в формуле для ). Эксперименты показывают (см. раздел 2), что для временного ряда длиной 100 периодов колебаний достаточно, чтобы абсолютная величина оценки АКФ не превышала примерно 0.2 при лагах, больших т.
4) Слабая индивидуальная нелинейность и слабый шум. Первое нужно для равномерности распределения разности фаз. Второе - для правомерности перехода от (6) к (7) с гауссовыми шумами . Эти условия можно сформулировать как
6£к ^ и \J~bk ^ . Во всех численных экспериментах в разделе 2 выполнялись
соотношения д£к < 0.2wk и < 0.1wk, чего оказывалось достаточно. Далее это условие подробнее не комментируется.
2. Численный эксперимент
Применимость предложенных оценок в случае временных рядов умеренной длины показана ниже в численных экспериментах для эталонных осцилляторов с различными свойствами. В п. 2.1 описана методика исследования, в пп. 2.2-2.5 -результаты расчетов для фазовых осцилляторов (п. 2.2), осцилляторов ван дер Поля (п. 2.3), систем Ресслера (п. 2.4) и систем Морриса-Лекара, являющихся моделью динамики нейрона (п. 2.5).
2.1. Методика исследования
В качестве эталонных примеров для тестирования предложенных оценок используются системы стохастических дифференциальных уравнений с различными видами нелинейности и силами связей. Уравнения численно интегрируются методом Эйлера с достаточно малым шагом h: h = 0.01 в пп. 2.2 и 2.3, h = 0.001 в п. 2.4, h = 0.1 в п. 2.5. Для корректного интегрирования стохастических дифференциальных уравнений нужно указать, в каком смысле понимается стохастический интеграл вида Jg(x(t))d%(t), где x - одна из переменных, ^ - белый шум, который входит в уравнения в виде слагаемого g(x) ■ ^ с гладкой функцией g. Допустимы различные варианты: интеграл Ито, интеграл Стратоновича, обобщенный стохастический интеграл, которые различаются при g(x) = const [40]. Однако во всех примерах, рассмотренных ниже, включая уравнения (6), используется g(x) = 1. При этом различные трактовки стохастического интеграла совпадают, так что согласно методу Эйлера в разностной схеме используется стохастическое слагаемое o-^Jh ■ en, где en - последовательность независимых гауссовых случайных величин с нулевым средним и единичной дисперсией.
Интервал выборки задается так, чтобы обеспечить примерно 10 точек на характерном периоде колебаний в пп. 2.2-2.4 (At = 0.6) и 25 точек на периоде в п. 2.5 (At = 1.0). Соответственно, величина т для оценки связей принимается равной т = 10At (пп. 2.2-2.4) и т = 25At (п. 2.5). Для каждой системы и каждого набора параметров генерируется большой ансамбль из Ns временных рядов, чтобы исследовать статистические свойства предложенных оценок связи. Во всех примерах принимается Ns = 100, если не указано другое значение. Начальные условия для получения различных временных рядов в ансамбле выбираются в соответствии
со стационарной плотностью распределения вектора состояния. Длина каждого временного ряда составляет примерно 100 характерных периодов, то есть N = 1000 в пп. 2.2-2.4, и N = 2500 в п. 2.5. Эту умеренную длину [7] уместно использовать как базовый случай. Зависимость результатов от N обсуждается кратко. По каждому ряду рассчитываются оценки связи и делаются выводы о том, обнаружено ли воздействие ] — к (положительный вывод) на заданном уровне значимости р или нет (отрицательный). Результаты расчетов представлены в основном для р = 0.05.
Обозначим Qj^k - число положительных выводов. Частота их появления
= Qj^k/Ns. Положительные выводы могут быть как «правильными» (если на самом деле есть воздействие ] — к), так и «ошибочными» (в противном случае). В каждом примере проверяется два обстоятельства. Во-первых, предложенные оценки можно считать применимыми только в тех случаях, когда частоты ошибочных положительных выводов не превышают заявленного уровня значимости: для всех пар (],к) таких, что Cj^k = 0, должно выполняться Vj^k(р) < Р с точностью до флуктуаций, определяющихся конечным размером ансамбля Ns. Эти флуктуации оцениваются следующим образом. Частота Vj^k - это оценка вероятности Pj^k вывода о наличии воздействия ] — к. По теореме Муавра-Лапласа величина Vj^k при большом Ns распределена по нормальному закону со средним Pj^k и дисперсией Pj^k(1 — Pj^k)/Ns. В частности, 95%-й интервал для Vj^k имеет вид Pj^k±1-96^/Pj^k(1 — Pj^k)/Ns. Поскольку вероятность ошибочных выводов Pj^k должна быть равна заявленному уровню значимости р, то Vj^k должна попадать в 95%-й интервал р ± 1.96^р(1 — p)/Ns. Применимость методики отвергается из-за высокой частоты ошибочных выводов, если оказывается, что Vj^k > р + 1-96\/p{Т—p)JN~s. Во-вторых, желательна максимальная чувствительность метода. Ее характеристика - частота правильных положительных выводов, то есть вероятность обнаружения существующих связей по отдельной временной реализации. Оценки сил воздействий и частот Vj^k рассчитываются для трех последовательно усложняющихся модельных функций (8): 1) «линейная» связь, то есть набор индексов 0.и, содержащий только т = п = 1; 2) допускается квадратичная нелинейность связи, то есть в 0.и есть пары индексов т = п = 1, т = 2, п = 1 и т = 1, п = 2; 3) допускается квадратичная или кубическая нелинейность связи, то есть в 0.и есть пары индексов т = п = 1, т = 2, п = 1, т = 1, п = 2, т = 3, п = 1 и т = 1, п = 3. Обозначим |0и| - размер множества 0.и. Тогда в первом случае |0и | = 1, то есть 0.и содержит одну пару индексов, во втором - |0и | =2, в третьем - |0и| = 3. Всегда принимается = 02. Число степеней свободы закона «хи-квадрат», используемого при оценке значимости выводов равно Ми = 2 |0и |.
2.2. Фазовые осцилляторы
Первый пример - фазовые осцилляторы (6) при различных видах и силах связи, соотношениях частот, интенсивностях шума.
2.2.1. Несвязанные фазовые осцилляторы. Случай несвязанных осцилляторов Си = 0 почти полностью соответствует теоретическим условиям применимости оценок. Однако они с гарантией применимы асимптотически, поэтому представляет интерес проверка практических результатов для рядов умеренной длины. Так, для
параметров о^ = 0.04, ю1 = 1.1, и ю2 = 0.9 результаты представлены на рис. 1. АКФ остаточных ошибок, оцененная по одному временному ряду (рис. 1, а), соответствует ожидаемой «треугольной» форме. Анализ ансамбля из N = 1000 рядов показывает, что величина р всегда не превышает 0.15, то есть весьма мала. На рис. 1, б показаны частоты положительных выводов о наличии связи для оценок с «линейной» моделью |Я&| = 1 в зависимости от уровня значимости р. Как и требуется, и У1^2 не превышают р. Те же результаты наблюдаются для моделей с |Я& | = 2 и |Я& | =3 (не показаны).
Трудности могут возникать при таких параметрах осцилляторов, когда оценки р часто принимают достаточно большие значения, то есть эмпирическое распределение разности фаз для отдельного временного ряда не является «почти равномерным». Это имеет место при уменьшении расстройки частот, что иллюстрирует рис. 2 для оценок с |Я& | = 1, где Ю1 = 1 + Дю, Ю2 = 1 — Дю. Так, при Дю = 0 значения р иногда превышают величину 0.5. На рис. 2, а показана доля временных рядов п0, для которых р < 0.45. При больших расстройках по = 1, а при малых уменьшается. На рис. 2, б представлены частоты положительных выводов без контроля условия р < рс для уровня значимости р = 0.05, а пунктиром показан допустимый уровень ошибок. При Дю частота ошибочных выводов для одного направления до-
Рис. 1. Несвязанные фазовые осцилляторы: а - АКФ остаточных ошибок £1 - оцененная по временному ряду (жирная линия) и теоретическая (тонкая линия); б - частоты положительных выводов о наличии воздействия 2 ^ 1 (сплошная линия) и 1 ^ 2 (пунктирная линия), полученные по ансамблю из 1000 временных рядов, обе частоты не превышают заданного уровня значимости р (штриховая линия)
Рис. 2. Несвязанные фазовые осцилляторы. Результаты оценивания с ¡Я*; | = 1 в зависимости от расстройки частот: а - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; б - частоты ошибочных положительных выводов 2 ^ 1 (жирная линия) и 1 ^ 2 (тонкая линия) и допустимый уровень ошибки (штриховая линия); в - то же, что на рис. б, но при дополнительном условии р < 0.45. Результаты оценивания с |Я*| =2 и |Я*| = 3 аналогичны (не показаны)
Рис. 3. Несвязанные фазовые осцилляторы. Результаты оценивания с |Яд;| = 1 в зависимости от уровня шума: а - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; б - частоты ошибочных положительных выводов 2 ^ 1 (жирная линия) и 1 ^ 2 (тонкая линия) и допустимый уровень ошибки (штриховая линия); в - то же, что на рис. б, при дополнительном условии р < 0.45 Результаты оценивания с |Яй| =2 и |Яй| = 3 аналогичны (не показаны)
стигает 0.12, что больше допустимого уровня. Если делать положительные выводы только при р < рс, то можно снизить вероятность ошибок, так как при достаточно малом рс лучше выполняется условие равномерного распределения разности фаз. В данном случае оказалось, что следует принять рс равным 0.45, чтобы обеспечить даже условную частоту ошибочных выводов = /(по^) не более допустимого уровня р + 1.96д/р(1 — р)/(п0М3) при любом Аю (рис. 2, в). Знаменатель по Ns - это количество временных рядов, удовлетворяющих условию р < рс. Условная частота ошибочных выводов не меньше полной частоты, так что требование ее малости строже, чем требование малости величины Qj^k/Ns.
Доля рядов с большими р растет при уменьшении и малой расстройке частот, например, при Аю = 0. Результаты представлены на рис. 3 для ^к| = 1, где вновь показано, что условие р < 0.45 дает условную частоту ошибочных положительных выводов не более допустимого уровня. Подбор рс был проведен и для всех примеров рассмотренных ниже. Значение рс = 0.45 оказалось достаточным во всех случаях, так что ниже все результаты приводятся только для оценивания с использованием критерия р < 0.45. Для ошибочных выводов ниже представлена их условная частота, а не полная.
2.2.2. Связанные фазовые осцилляторы. Рассмотрим связанные фазовые осцилляторы (6) с Ю1 = 1.1, Ю2 = 0.9, а^ = 0.04 и однонаправленной связью 2 ^ 1, то есть С2 = 0. Начнем с «первого порядка» нелинейности: С1(ф2, Ф1) = = ^2^1 йш(ф2 — Ф1). Положительные выводы о наличии влияния 1 ^ 2 ошибочны и для применимости методики их вероятность должна быть не более р. Выводы о наличии связи 2 ^ 1 при ^2^1 > 0 правильные, так что желательна их максимальная вероятность. Условие на АКФ остаточных ошибок еи выполняется, так как она аналогична рис. 1, а (не показана).
Рис. 4 иллюстрирует результаты оценивания при увеличении коэффициента связи к2^1. Доля рядов с р < 0.45 уменьшается с ростом к2^1 (рис. 4, а). Однако частота ошибочных выводов остается не более допустимой при любом к2^1 (рис. 4, б). Это имеет место для любого из трех различных ^к |. Частота правильных выводов здесь и во всех примерах ниже рассчитывается в отношении ко всему размеру ансамбля: у2^1 = Q2^1/Ns без деления на по, то есть это оценка полной вероятности обнаружения воздействия (не условной). При малых к2^1 чувствитель-
По-.
0.8
0.4 :
0: 0
а
Рис. 4. Оценка однонаправленной связи между фазовыми осцилляторами в случае функции связи С].(ф2, Ф1) = к2^1 вш(ф2 — Ф1): а - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; б - частоты ошибочных выводов о наличии воздействия 1 ^ 2 (жирная линия - для = 1, тонкая - для | = 2, штриховая - для | = 3, пунктир - допустимый уровень ошибки); в - частоты правильных выводов о наличии воздействия 2 ^ 1 (те же обозначения, что и на рис. б; г - средние значения оценок сил воздействия 2 ^ 1 (кружки, соединенные линией), и 1 ^ 2 (крестики)
ность метода растет с ростом (рис. 4, в). При этом чувствительность выше для меньшего |. Это ожидаемый результат, так как связь между осцилляторами в данном примере описывается многочленом вида 8т(ф1 — ф2) и отражается уже моделью с | = 1. Увеличение | приводит только к появлению «лишних» оцениваемых коэффициентов, то есть к росту дисперсии оценок без улучшения описания динамики исходной системы. Это и является причиной снижения процента статистически значимых выводов о наличии воздействия. Использование критерия р < 0.45 ведет к отсеиванию некоторых временных рядов ансамбля и снижает чувствительность метода при больших Средние по ансамблю значения оценок С^^к показаны на рис. 4, г. Сила существующего воздействия (р2^1 растет квадратично с ростом к2^1 (кружки). Это ожидаемый результат, так как по определению (10) эта характеристика есть сумма квадратов коэффициентов связи. Оценка (р1^2 в среднем равна нулю (крестики). Величины (с!^^^ одинаковы для всех | и их графики полностью накладываются друг на друга на рис. 4, г.
На рис. 5 и 6 аналогично представлены результаты для связи, заданной нерезонансными членами С1(ф2, ф1) = к2^1 8ш(2ф2 — ф1) и С1(ф2, ф1) = к2^1 8т(3ф2 — Ф1), соответственно. Это связи более высокого порядка нелинейности, прочие параметры осцилляторов те же. Отличие этих случаев в том, что р даже не превышает величины 0.15, так что можно использовать все временные ряды ансамбля для оценивания (рис. 5, а и 6, а). Частота ошибочных выводов не превышает допустимой при использовании любого | (рис. 5, б и 6, б). Вероятность правильных выводов для связи 8ш(2ф2 — Ф1) равна нулю при использовании модели с низким порядком нелинейности | = 1 (рис. 5, в). Для | =2 и = 3 чувствительность метода растет с ростом к2^1. При этом для | =3 она ниже, так как модель фазовой динамики содержит лишние члены по сравнению с | = 2. Величина С 1^2 в среднем равна нулю для всех моделей (не показана). Среднее значение С2^1 представлено на рис. 5, г. Для | = 1 эта величина равна нулю (жирная линия), так как модель не может описать связь более высокого порядка нелинейности. Для двух других моделей одинакова и демонстрирует ожидаемую квадратичную зависимость от к2^1 (тонкая линия). Аналогично, связь эт(3ф2 — Ф1) обнаруживается только при использовании | = 3, а более простые модели не могут ее «почувствовать» (рис. 6, в, г).
Рис. 5. Оценка однонаправленной связи между фазовыми осцилляторами в случае функции связи С].(ф2, Ф1) = к2^1 вш(2ф2 — Ф1): а - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; б - частоты ошибочных выводов о наличии воздействия 1 ^ 2 (жирная линия - для | = 1, тонкая - для |Як| = 2, штриховая - для | = 3, пунктир - допустимый уровень ошибки); в - частоты правильных выводов о наличии воздействия 2 ^ 1 (те же обозначения, что и на рис. б; г - среднее по ансамблю значение оценки силы воздействия 2 ^ 1 для модели с | = 1 (жирная линия) и | = 2 (тонкая линия). Для |Як| =3 график в точности такой же, как для |Як| =2
Рис. 6. Оценка однонаправленной связи между фазовыми осцилляторами в случае функции связи С1(ф2, ф1) = к2^1 вш(3ф2 — ф1): а - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; б - частоты ошибочных положительных выводов о наличии воздействия 1 ^ 2 (жирная линия - для | = 1, тонкая - для | = 2, штриховая - для |Як| = 3, пунктир - допустимый уровень ошибки); в - частоты правильных положительных выводов о наличии воздействия 2 ^ 1 (те же обозначения, что на рис. б; г - среднее по ансамблю значение оценки силы воздействия 2 ^ 1 для модели с | = 1 (жирная линия) и = 3 (тонкая линия). Для =2 график такой же, как для = 1
Аналогичные результаты получены при соотношении частот примерно равном 1:2 и 1:3 и различных порядках нелинейности связи, при однонаправленном воздействии в любую сторону и при двунаправленной связи. Во всех случаях предложенные оценки применимы. При этом результаты сохраняются при уменьшении длины ряда до 50 характерных периодов, и только еще меньшая длина ряда приводит к большему числу ошибочных выводов из-за недостаточной статистики (графики не показаны).
Заметим, что в случае двунаправленной связи трудности оценивания меньше, поскольку ошибочных положительных выводов не существует из-за наличия воздействий в обе стороны, а опасность есть только в получении смещенных оценок сил воздействия. Однако при достаточно слабых связях (что контролируется по величине р) оценки Cj^k оказываются несмещенными. Графики аналогичны полученным в статье [7, рис. 3] и здесь не приводятся. То же наблюдение имеет место и для других осцилляторов. Поэтому ниже иллюстрируется только более сложный случай однонаправленной связи, где нужно контролировать частоту ошибочных выводов о наличии воздействий.
2.3. Осцилляторы ван дер Поля
Обратимся к примерам, где ситуация усложняется из-за необходимости расчета фаз по временным рядам и изменения свойств осцилляторов, в частности, свойств шумов ей в уравнениях фазовой динамики. Так, однонаправленно связанные осцилляторы ван дер Поля задавались уравнениями
^-Х- — (0.2 — х\) + ®2Х1 = + к2^1д(х2),
й2х2
— (0.2 — х2) ^ + «2x2 =
(15)
где ^1)2 - независимые источники белого шума, д(х2) - функция связи.
В отсутствие шумов осцилляторы демонстрируют автоколебания с круговыми частотами «1,2, и для каждого отдельного осциллятора справедливо приближение медленно меняющихся (по сравнению с периодом колебаний) амплитуд. При этом фазовая динамика описывается уравнениями фазовых осцилляторов (6) с Ок = 0. Наличие связи отражается в появлении дополнительных членов в уравнении фазовой динамики для первого осциллятора. Если д(х2) = х2, то получим О1(ф1, ф2) а: к2^1 8ш(ф2 — ф1), учитывая, что при отсутствии шума х2 а ехр(гф2(£)). Если д(х2) = х2, то есть д а ехр(йф2(£)), получим связь более высокого порядка нелинейности: О1(ф1, ф2) а к2^1 8т(2ф2 — ф1). Аналогично, при д(х2) = х] — 3х2/4 имеем д а ехр(г3ф2(£)) и О1(ф1, ф2) а к2^1 эш(3ф2 — ф1).
При сильных шумах ^12 амплитуды колебаний осцилляторов заметно флуктуируют, то есть не являются медленно меняющимися. Поэтому все рассуждения применимы только приближенно. Такой случай и рассмотрим в качестве теста для предложенных оценок связи. Параметры осцилляторов 0| = 0.2, « = 1.05, коэффициент связи меняется от 0 до 0.2, связь резонансная, что обеспечивается различными порядками нелинейности для различных Ю2:
1) ю2 = 0.95 (соотношение частот ю1 : ю2 примерно 1:1) и д(х2) = х2,
2) ю2 = 0.5 (соотношение частот примерно 2:1) и д(х2) = х2,
3) ю2 = 0.35 (соотношение частот примерно 3:1) и д(х2) = х] — 3х2/4. Фазы рассчитывались по сигналам х1>2(£) с помощью преобразования Гильберта. Рис. 7, а показывает, что, несмотря на заметные флуктуации амплитуды, чередование максимумов и минимумов сигнала с основным периодом четко видно, так что фазы хорошо определены. АКФ остаточных ошибок ек имеет требуемую форму (рис. 7, б).
Рис. 7. Осцилляторы ван дер Поля (15): а - фрагмент временной реализации первого осциллятора при отсутствии связи; б - АКФ остаточных ошибок е1 - оценка (жирная линия) и требуемая теоретически (тонкая линия)
Рис. 8. Оценка однонаправленной связи между осцилляторами ван дер Поля (15). Первая строка - для функции связи д(х2) = х2, вторая - для д(х2) = х2, третья - для д(х2) = х2 — 3х2/4. Первый столбец - доля временных рядов с рт,п < 0.45 в ансамбле из 100 рядов. Второй - частоты ошибочных выводов о наличии воздействия 1 ^ 2 (жирная линия - для = 1, тонкая - для | = 2, штриховая - для | = 3, пунктир - допустимый уровень ошибки). Третий - частоты правильных выводов о наличии воздействия 2 ^ 1. Четвертый - среднее по ансамблю значение оценки силы воздействия 2 ^ 1 для модели с | = 1 (жирная линия), | = 2 (тонкая линия), и | = 3 (штриховая линия)
Результаты оценивания, демонстрирующие применимость методики, представлены на рис. 8. Первый столбец показывает долю рядов с р < 0.45 для линейной связи (рис. 8, а), р1>2 < 0.45 для квадратичной (рис. 8, д) и р1)3 < 0.45 для кубической (рис. 8, и). Всегда частота ошибочных выводов не более допустимой (рис. 8, второй столбец) для любого |Ои |. При этом во втором случае (рис. 8, е) на длине временного ряда укладывается только 50 характерных периодов второго осциллятора, а в третьем (рис. 8, к) - 30 периодов. Тем не менее, методика не дает большего числа ошибок, что иллюстрирует ее применимость даже для таких коротких рядов. Вероятность обнаружения существующей связи растет с ростом к2^1 (рис. 8, третий столбец). В первом случае достаточно |Ои | = 1, а большие значения лишь снижают чувствительность (рис. 8, в, г). Во втором случае лучше всего [□к| = 2, тогда как |Ои| = 1 не обнаруживает квадратичную связь (рис. 8, ж, з). В третьем случае кубическая связь приводит к тому, что в уравнениях фазовой динамики есть и линейные и кубические члены. Поэтому методика для |Ои | = 1 и |Йк | = 3 примерно одинаково чувствительна (рис. 8, л). Причина состоит в том, что кубическая модель лучше описывает нелинейность (это показывают и оценки силы связи на рис. 8, м), но содержит больше оцениваемых коэффициентов. Сов-
местное влияние этих двух факторов приводит к той же чувствительности, что и для [□к| = 1. Использование |Ок| =2 в основном снижает чувствительность, так как квадратичные члены в уравнении фазовой динамики отсутствуют или очень малы (см. рис. 8, м).
2.4. Системы Ресслера
Системы Ресслера также рассматривались с однонаправленной связью
йх1/(Л, = — «1у1 — г1,
йу1/(И = «х + ау1 + к2^1д(у2) + (16)
(г1/М = Ь + (х1 — е)г1,
(х2/(Ь = —«2У2 — ¿2,
((у2/(И = «2х1 + ау2 + ^2, (17)
(г2/(Ь = Ь + (х2 — с)г2,
где ^1,2 - независимые источники белого шума, д(у2) - функция связи. Выбирались такие значения а, Ь, с, при которых осцилляторы (при 0| = 0 и к2^1 = 0) демонстрируют колебания с основными частотами «1,2 в режиме спирального хаоса. Результаты анализа этих систем оказываются во многом схожими с предыдущими примерами. Поэтому далее представлены иллюстрации для «1 = 1.05, «2 = 0.95, д(у2) = у2 и следующих ситуаций:
1) а = 0.2, Ь = 0.2, с = 6.1 [41] и 0| = 0.1,
2) а = 0.398, Ь = 2.0, с = 4.0 [42] и 0| = 0.2, д(у2) = у2.
При отсутствии шумов ^12 в (16) фазовую динамику систем Ресслера можно описать уравнениями фазовых осцилляторов (6), где шумы ^12 в (6) присутствуют как результат влияния хаотически меняющихся амплитуд. Эти шумы не белые, что дает примеру новизну по сравнению с предыдущими.
При анализе рассматривались сигналы х1;2. Фазы хорошо определены (рис. 9, а и 10, а) и рассчитывались с помощью преобразования Гильберта. АКФ остаточных ошибок ек для первого набора параметров имеет почти требуемую форму, хотя и не совсем: АКФ принимает значение 0.4 при некоторых лагах, больших т (рис. 9, б). Однако такое отклонение оказывается допустимым, так как результаты оценивания вполне аналогичны представленным выше (рис. 9, в-е): методика не дает ошибочных выводов сверх допустимой величины (см. рис. 9, г); слишком большое значение |Йк | приводит к снижению чувствительности (см. рис. 9, д, е).
Для второго набора параметров несмотря на хорошо определенную фазу, АКФ остаточных ошибок ек сильно отличается от требуемой «треугольной» функции (рис. 10, б). Величина АКФ достигает 0.8 для многих лагов, значительно больших т. Как результат, число ошибочных выводов о наличии воздействия 1 ^ 2 гораздо больше допустимого (рис. 10, г), то есть методика не применима. Таким образом, выполнение условия на АКФ остаточных ошибок оказывается весьма существенным.
Связь по переменной y в данном примере - это аналог связи по dx/dt для осцилляторов ван дер Поля, что можно увидеть, если записать уравнение для d2x/dt2, исключив y из первого уравнения системы Ресслера. Это также отличает данный пример от предыдущего. Связь по x проверялась таким же образом, результаты оказались аналогичны (не показаны).
В работе выборочно проверялись возможности использования не «одинаковых» переменных для введения фаз. Например, xi для первой системы и y2 или Z2 (вместо x2) - для второй. Результаты оказываются вполне аналогичными. При использовании y2 значения фазы ф2, рассчитанные с помощью преобразования Гильберта, только сдвигаются на п/2 по сравнению с фазами, рассчитанными по сигналу x2, что не влияет на характеристики связи. Координата z2 демонстрирует поведение импульсного характера, поэтому для расчета фазы предпочтительнее пороговый метод (см. в п. 2.5). Он дает фазу, которая является нелинейной функцией величины ф2, полученной по сигналу x2. Но эта нелинейная функция очень близка к линейной H(Ф2) ~ Ф2 + const, то есть значения фазы слабо меняются и все результаты оценивания остаются теми же самыми (не показаны).
2.5. Системы Мориса-Лекара
Наконец, рассмотрим пример систем с «импульсным» характером воздействия. Это системы Мориса-Лекара [43], моделирующие динамику нейрона,
d" = h - gi(xi - Vi) - gKVi(xi - Vk) - gcam^(xi - Vea) + k^i/(xi,x2), d- = X(xi) (w^(xi) - yi) + h,
at (18)
= I2 - gi(x2 - Vi) - дкУ2(x2 - Vk) - gcam^(x2 - Vca),
d2 = ^(x2) (w^(x2) - У2) + Ъ,
где m^(x) = (1+tanh((x - Vi)/V2))/2, w^(x) = (1 + tanh((x - Vs)/V4))/2, X(x) = (cosh ((x - Vs)/2V4))/3, Vi = -0.1, V2 = 0.15, V3 = 0.1, V4 = 0.145, Vi = -0.5, Vk = -0.7, Vea = 1.0, gi = 0.5, gK = 2.0, gca = 1.33, Ii = I2 = 0.075. При отсутствии шума и связи системы демонстрируют периодические автоколебания [44,45]. Введение слабого шума с 0| = 0.005 слабо возмущает эти периодические колебания [46]. Связь пороговая, что имитирует свойства взаимодействия
нейронов: /(x-,x2) = (1 - xi) ( -1 + —-г1-—— ) , x2 > 0, /(x-,x2) = 0,
2 1+exp(-Ж2/0.1)< Х2 < 0. Фазовая динамика таких «нейроноподобных» систем может приближенно описываться уравнениями фазовых осцилляторов (6), но свойства шумов в них аналитически получить не удается, а связь может быть существенно нелинейной и задаваться даже дельта-функцией [45], поскольку импульсный характер обнаруживают временные реализации (рис. 11, а).
Рис. 9. Системы Ресслера (16), первый набор параметров: а - фрагмент временной реализации первого осциллятора при отсутствии связи; б - АКФ остаточных ошибок е1 - оценка (жирная линия) и требуемая теоретически (тонкая линия); в - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; г -частоты ошибочных выводов о наличии воздействия 1 ^ 2 (жирная линия - для \ = 1, тонкая - для =2, штриховая - для \ = 3, пунктир - допустимый уровень ошибки); д - частоты правильных выводов о наличии воздействия 2 ^ 1; е - среднее по ансамблю значение оценки силы воздействия 2 ^ 1 для | = 1 (жирная линия), | = 2 (тонкая линия), и = 3 (штриховая линия сливается с тонкой)
Рис. 10. Системы Ресслера (16), второй набор параметров, те же обозначения, что на рис. 9: а -фрагмент временной реализации; б - АКФ остаточных ошибок е1; в - доля рядов с р < 0.45 в ансамбле из 100 рядов; г - частоты ошибочных выводов о наличии воздействия 1 ^ 2; д - частоты правильных положительных выводов о наличии воздействия 2 ^ 1; е - среднее по ансамблю значение оценки силы воздействия 2 ^ 1 (для всех одинаково)
О 0.004 0.008 г 0 0.004 0.008 ¿) 0 0.004 0.008 е 0 0.004 0.008
Рис. 11. Системы Мориса-Лекара (17), те же обозначения, что на рис. 9, 10: а - фрагмент временной реализации, пунктиром показан уровень хх = 0 для введения фазы; б - АКФ остаточных ошибок £х; в - доля временных рядов с р < 0.45 в ансамбле из 100 рядов; г - частоты ошибочных выводов о наличии воздействия 1 ^ 2; д - частоты правильных выводов о наличии воздействия 2 ^ 1; е -среднее по ансамблю значение оценки силы воздействия 2 ^ 1 (для =2 и | = 3 совпадают)
При анализе рассматривались сигналы х\,2. Фазы рассчитывались с помощью линейной интерполяции: фаза сигнала линейно нарастает на 2п между двумя последовательными пересечениями реализацией Хк(¿) оси абсцисс (см. рис. 11, а, пунктирная линия) снизу вверх. Фазы хорошо определены [46]. Оказывается, что АКФ остаточных ошибок Ек имеет требуемую форму (рис. 11, б), так что это теоретическое условие применимости методики выполняется.
Результаты оценивания связи (рис. 11, в-е) аналогичны представленным выше примерам. Несмотря на сложный вид функции связи в исходных уравнениях (19) и ее неизвестную форму в модели (6) или (7), наличие и интенсивность связи адекватно определяются.
Заключение
В работе предложены новые характеристики направленной связи между осцилляторами и их оценки по временным рядам, основанные на идее моделирования фазовой динамики. В отличие от ранее использовавшихся, предложенные характеристики пригодны в случае произвольного порядка нелинейности связей. Аналитически выведено выражение для уровня статистической значимости, на котором делается вывод о наличии воздействий. Это важно для приложений, так как обеспечивает заданную надежность при работе с относительно короткими зашумленными временными рядами.
Представлены количественные условия применимости методики на практике: хорошо определенные фазы колебаний, невысокая степень синхронности сигналов, ограничения на корреляционные свойства и интенсивность шумов в уравнениях фазовой динамики. Работоспособность методики показана в численных экспериментах на эталонных осцилляторах. Показано, что для описания более высокого порядка
нелинейности требуется использовать более сложную модель фазовой динамики. Однако использование сложной модели без необходимости только снижает чувствительность метода.
Строго говоря, развитый формализм применим асимптотически, но эксперименты показывают, что практические приложения возможны для рядов умеренной длины. В работе рассматривались в основном ряды длиной 100 характерных периодов, этого оказалось достаточно во всех случаях. Более того, для рассмотренных простых моделей фазовой динамики (описывающих линейные, квадратичные или кубические члены воздействия) длина ряда может быть уменьшена до 30-50 характерных периодов, что иллюстрируют некоторые представленные примеры. Однако для еще более сложных моделей, содержащих много членов более высоких порядков, минимальная требуемая длина ряда должна быть больше и растет с ростом сложности модели. Детальное количественное исследование закона этого роста заняло бы слишком много времени. Кроме того, следует ожидать, что на практике чаще нужно извлекать информацию из небольших объемов данных, для чего придется ограничиваться использованием невысоких порядков нелинейности в моделях. Если все же требуется использовать модель очень высокого порядка, то необходимо провести предварительные численные тесты и установить требования на минимальную длину ряда.
В заключение отметим, что эффективность методики показана для осцилляторов с различными свойствами (фазовые осцилляторы, осцилляторы ван дер Поля, системы Ресслера, системы Мориса-Лекара), то есть различными функциями связи и свойствами шумов в уравнениях фазовой динамики. Это свидетельствует о ее значительной степени общности и возможности широкого практического использования.
Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант 08-02-00081), программ РАН и Министерства образования и науки РФ.
Библиографический список
1. Пиковский А.С., Розеблюм М.Г., Куртс Ю. Синхронизация: фундаментальное нелинейное явление. М.: Техносфера. 2003.
2. Tass P.A. Phase resetting in medicine and biology - stochastic modelling and data analysis. Berlin: Springer, 1999.
3. Boccaletti S., Kurths J., Osipov G., Valladares D., Zhou C. The synchronization of chaotic systems // Phys. Rep. 2002. Vol. 366. P. 1.
4. Mosekilde E., Maistrenko Yu., Postnov D. Chaotic Synchronization. Applications to Living Systems. Singapore: World Scientific, 2002.
5. Balanov A., Janson N., Postnov D., Sosnovtseva O. Synchronization: From Simple to Complex. Berlin: Springer-Verlag, 2008.
6. Rosenblum M.G., Pikovsky A.S. Detecting direction of coupling in interacting oscillators // Phys. Rev. E. 2001. Vol. 64, 045202(R).
7. Smirnov D.A., Bezruchko B.P. Estimation of interaction strength and direction from short and noisy time series // Phys. Rev. E. 2003. Vol. 68, 046209.
8. Palus M., Stefanovska A. Direction of coupling from phases of interacting oscillators: An information-theoretic approach // Phys. Rev. E. 2003. Vol. 67, 055201.
9. Kralemann B., Cimponeriu L., Rosenblum M., Pikovsky A., Mrowka R. Uncovering interaction of coupled oscillators from data// Phys. Rev. E. 2007. Vol. 76, 055201.
10. Kralemann B., Cimponeriu L., Rosenblum M., Pikovsky A., Mrowka R. Phase dynamics of coupled oscillators reconstructed from data // Phys. Rev. E. 2008. Vol. 77, 066205.
11. Baccala L.A., Sameshima K. Partial directed coherence: a new concept in neural structure determination // Biological Cybernetics, 2001. Vol. 84. P. 463.
12. Ancona N., Marinazzo D., Stramaglia S. Radial basis function approach to nonlinear Granger causality of time series // Phys. Rev. E. 2004. Vol. 70, 056221.
13. Schreiber T. Measuring information transfer // Phys. Rev. Lett. 2000. Vol. 85. P. 461.
14. Verdes P.F. Assessing causality from multivariate time series // Phys. Rev. E. 2005. Vol. 72, 026222.
15. Hlavackova-Schindler K., Palus M., Vejmelka M., Bhattacharya J. Causality detection based on information-theoretic approaches in time series analysis // Phys. Rep. 2007. Vol. 441. P. 1.
16. Vejmelka M., Palus M. Inferring the directionality of coupling with conditional mutual information // Phys. Rev. E. 2008. Vol. 77, 026214.
17. Rosenblum M.G., Cimponeriu L., Bezerianos A., PatzakA., Mrowka R. Identification of coupling direction: Application to cardiorespiratory interaction // Phys. Rev. E. 2002. Vol. 65, 041909.
18. Prokhorov M.D., Ponomarenko V.I., Gridnev V.I., Bodrov M.B., Bespyatov A.B. Synchronization between main rhythmic processes in the human cardiovascular system // Phys. Rev. E. 2003. Vol. 68, 041913.
19. Luchinsky D.G., Millonas M.M., Smelyanskiy V.N., Pershakova A., Stefanovska A., McClintock P. V. Nonlinear statistical modeling and model discovery for cardiorespiratory data // Phys. Rev. E. 2005. Vol. 72, 021905.
20. Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. Detection of synchronization from univariate data using wavelet transform // Phys. Rev. E. 2007. Vol. 75, 056207.
21. Sosnovtseva O.V., Pavlolv A.N., Mosekilde E., Holstein-Rathlou N.H., Marsh D.J. Double-wavelet approach to study frequency and amplitude modulation in renal autoregulation // Phys. Rev. E. 2004. Vol. 70, 031915.
22. Pavlov A.N., Sosnovtseva O.V., Pavlova O.N., Mosekilde E., Holstein-RathlouN.-H. Characterizing multimode interaction in renal autoregulation // Physiological Measurements, 2008. Vol. 29. P. 945.
23. Eguia M.C., Dawson S.P., Mindlin G.B. Information transmission and recovery in neural communications channels // Phys. Rev. E. 2000. Vol.62. P. 7111.
24. Kiemel T., Gormley K., Guan L., Williams T., Cohen A. Estimating the strength and direction of functional coupling in the lamprey spinal cord // J. Computational Neuroscience. 2003. Vol. 15. P. 233.
25. Blinowska K.J., Kus R., Kaminski M. Granger causality and information flow in multivariate processes // Phys. Rev. E. 2004. Vol. 70, 050902(R).
26. Pereda E., Quian Quiroga R., Bhattacharya J. Nonlinear multivariate analysis of neurophysiological signals // Progress in Neurobiology. 2005. Vol. 77. P. 1.
27. Brea J., Russell D.F., Neiman A.B. Measuring direction in the coupling of biological oscillators: A case study for electroreceptors of paddlefish // Chaos. 2006. Vol. 16, 026111.
28. Schelter B., Winterhalder M., Eichler M., Peifer M., Hellwig B., Guschlbauer B., Luecking C., Dahlhaus R., Timmer J.Testing for directed influences among neural signals using partial directed coherence // J. Neurosci. Methods. 2005. Vol. 152. P. 210.
29. Wang S., Chen Y., Ding M., Feng J., Stein J.F., Aziz T.Z., Liu XJ.Revealing the dynamic causal interdependence between neural and muscular signals in Parkinsonian tremor // J. Franklin Institute. 2007. Vol. 344. P. 180.
30. Osterhage H., Mormann F., Wagner T., Lehnertz K. Measuring the directionality of coupling: Phase versus state space dynamics and application to EEG time series // Int. J. Neural Syst. 2007. Vol. 17. P. 139.
31. Smirnov D.A., Barnikol U.B., Barnikol T.T., Bezruchko B.P., Hauptmann C., Buehrle C., Maarouf M., Sturm V., Freund H.-J., Tass P.A. The generation of Parkinsonian tremor as revealed by directional coupling analysis // Europhys. Lett. 2008. Vol. 83, 20003.
32. Mokhov I.I., Smirnov D.A. El Nino Southern Oscillation drives North Atlantic Oscillation as revealed with nonlinear techniques from climatic indices // Geophysical Research Letters. 2006. Vol. 33, 024557.
33. Smirnov D.A., Bezruchko B.P. Detection of couplings in ensembles of stochastic oscillators // Phys. Rev. E. 2009. Vol. 79, 046204.
34. Kori H., Kuramoto Y Slow switching in globally coupled oscillators: robustness and occurrence through delayed coupling // Phys. Rev. E. 2001. Vol. 63, 046214.
35. Kuramoto Y Chemical Oscillations, Waves and Turbulence. Berlin: Springer-Verlag, 1984.
36. Pikovsky A.S., Rosenblum M.G., Kurths J.Phase synchronization in regular and chaotic systems // Int. J. Bifurc. Chaos. 2000. Vol.10, № 10. P. 2291.
37. Gabor D. Theory of communication // J. Inst. Elect. Eng. (London). 1946. Vol. 93. P. 429.
38. Rosenblum M.G., Pikovsky A.S., Kurths J., Schaefer C., Tass P.A. Phase synchronization: from theory to data analysis // Neuro-informatics. Handbook of Biologi cal Physics. Edited by F. Moss and S. Gielen. New York: Elsevier Science, 2001. P. 279.
39. Kendall M.C., Stuart A. The advanced theory of statistics. New York: Hafner, 1979.
40. Никитин Н.Н., Разевиг В.Д. Методы цифрового моделирования стохастических дифференциальных уравнений и оценка их погрешностей // Журнал выч. математики и мат. физики. 1978. Т. 18, № 1. С. 106.
41. Анищенко В.С., Вадивасова Т.Е., Окрокверцхов Г.А., Стрелкова Г.И. Статистические свойства динамического хаоса // УФН. 2005. Т. 175. С. 163.
42. Thompson J.M.T. and Stewart H.B. Nonlinear Dynamics and Chaos. New York: Wiley, 1987.
43. Morris C., LecarH. Voltage oscillations in the barnacle giant muscle fiber // Biophys. J. 1981. Vol. 35. P. 193.
44. Izhikevich E.M. Neural excitability, spiking and bursting // Int. J. Bifurc. Chaos. 2000. Vol. 10. P. 1171.
45. Ermentrout G.B., Kopell N. Oscillator death in systems of coupled neural oscillators // SIAM J. Appl. Math. 1990. Vol. 50. P. 125.
46. SmirnovD., SchelterB., Winterhalder M., TimmerJ. Revealing direction of coupling between neuronal oscillators from time series: Phase dynamics modeling versus partial directed coherence // Chaos. 2007. Vol. 17, 013111.
Саратовский филиал ИРЭ Поступила в редакцию 29.09.2009
им. В.А. Котельникова РАН После доработки 3.03.2010
REVEALING NONLINEAR COUPLINGS BETWEEN STOCHASTIC OSCILLATORS FROM TIME SERIES
D.A. Smirnov
The problem of detection and quantitative characterization of nonlinear directional couplings between stochastic oscillators is considered. Coupling characteristics and a technique for their estimation from time series are suggested. An analytic expression for a statistical significance level of the conclusion about coupling presence is derived that allows a reliable inference from relatively short signals. Performance of the approach is demonstrated in numerical experiments with diverse individual properties of oscillators and different kinds of coupling functions.
Keywords: Coupling estimation, coupled oscillators, nonlinear time series analysis, phase dynamics, statistical inference
Смирнов Дмитрий Алексеевич - родился в 1977 году, окончил факультет нелинейных процессов Саратовского госуниверситета (1999), защитил кандидатскую диссертацию (2001). Старший научный сотрудник СФИРЭ РАН. Опубликовал более 40 статей в научных журналах и 2 монографии (в соавторстве). Область научных интересов: теория колебаний и волн, теория динамических систем, анализ временных рядов, математическое моделирование сложных систем по данным наблюдений. 410019, Саратов, ул. Зеленая, д. 38 СФИРЭ им. В.А. Котельникова РАН E-mail: [email protected], [email protected]
