CC BY
44
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
Изв. вузов «ПНД», т. 21, № 2, 2013 УДК 530.182,51 -73
ВЛИЯНИЕ ВЫБОРА СТРУКТУРЫ МОДЕЛИ НА РАБОТОСПОСОБНОСТЬ МЕТОДА НЕЛИНЕЙНОЙ ПРИЧИННОСТИ ПО ГРЕЙНДЖЕРУ
М. В. Корнилов, И. В. Сысоев
В настоящее время метод нелинейной причинности по Грейнджеру активно используется в ряде приложений медицины, биологии, физики для выявления направленной связанности между объектами по записям их колебаний (временным рядам) с помощью предсказательных моделей. В работе исследуется влияние выбора структуры модели на работоспособность метода. Численно на примере связанных эталонных уравнений продемонстрирована возможность получения достоверных оценок, даже если структура предсказательной модели отличается от структуры эталонной системы.
Ключевые слова: Метод нелинейной причинности по Грейнджеру, реконструкция по временным рядам, нелинейные динамические системы.
Введение
Причинность по Грейнджеру - известный способ оценки влияния одной системы на другую по записям их колебаний (временным рядам), впервые предложенный нобелевским лауреатом по экономике К. Грейнджером [1]. Согласно данному методу степень влияния одной системы на другую оценивается по изменению точности прогноза поведения первой системы при введении в прогностическую математическую модель данных о колебаниях второй системы. Уменьшение ошибки прогноза истолковывается как признак влияния второй системы на первую. Метод неоднократно применялся для оценки связанности экспериментальных систем, например, между различными отведениями энцефалограммы головного мозга [2-4], процессами в климатологии [5,6]. В указанных случаях результаты получались с помощью авторегрессионных моделей с полиномиальными [7, 8] или радиальными [9] базисными функциями, выбор которых опирался на физиологические и физические соображения, а их достоверность подтверждалась анализом характеристик остатков и с помощью суррогатных данных.
Но не случаен ли успех? Остаётся открытым вопрос влияния вида предсказательной модели на оценку причинно-следственных связей методом причинности
по Грейнджеру. Ответ на него целесообразно искать, используя в качестве объектов эталонные связанные уравнения, так как модели реальных объектов, как правило, не единственны и вопрос их адекватности требует доказательств. При таком выборе объекта нам известна его структура, включая наличие и характер связей, что позволяет проверить правильность оценок с помощью метода причинности по Грейнджеру.
Целью данной работы является проверка на тестовых примерах (связанных эталонных динамических системах) работоспособности метода при различных уровнях связи и значимости результатов причинности по Грейнджеру в случае, когда структура модели отлична от структуры объекта. Рассматривается модификация метода с полиномиальными нелинейными функциями.
1. Описание метода
Пусть имеются два временных ряда: ряд |жга}^=1 от системы X и ряд {уп}^= 1 от системы У, где п - дискретное время, N - длина временного ряда. На основе анализа реализаций {хп}^=1 и {уп}^=1, которые в общем случае содержат и шумы, требуется определить, влияет ли система У на систему X или нет.
На первом шаге строится индивидуальная модель
хП = /(Хп-1,Хп-21, Хп-ВэI, С), (1)
где / - аппроксимирующая функция, I - лаг модели, Ds - собственная размерность модели, ^ - неизвестный вектор коэффициентов. В качестве функции / можно использовать разложение по некоторому базису, например, степенной многочлен общего вида. В основополагающей работе [1] используются линейные функции. Коэффициенты ^ подбираются методом наименьших квадратов по временному ряду {хп}п=1. Среднеквадратичная ошибка аппроксимации построенной модели определяется по формуле = ^ (х'п — хп)2.
На втором шаге строится совместная модель
хп = д(Хп-1 ,Хп-21, ...,Хп--ОвиУп-1-т, ...,Уп-Ва1-т, С ), (2)
где для оценки коэффициентов аппроксимирующей функции д используется также реализация системы У; Da размерность добавки (число учитываемых значений из ряда {уп}^=1); С - коэффициенты совместной модели; т - возможная задержка, обусловленная конечным временем распространения сигнала (на практике часто неизвестна и должна быть определена).
Построив совместную модель, рассчитывают среднеквадратичную ошибку прогноза совместной модели е2- = ^ {х"п — .г,, )2. При е2- < говорят, что У действует на X (системы связаны).
В качестве меры связанности используется показатель улучшения прогноза
Р1 = 1-4-. (3)
4
Если Р1 = 0 (учёт сигнала У не помог в предсказании X), считают, что У не воздействует на X. Если же Р1 ^ 1 (учёт сигнала У существенно улучшил предсказание X), считают, что У воздействует на X.
В реальной ситуации всегда присутствуют шумы, длина временного ряда ограничена, точность вычислений конечна и т.п., поэтому требуется проводить проверку значимости выводов - оценку вероятности, что полученный результат не случаен. Известны различные виды такого тестирования: теоретические оценки [10], использование суррогатных данных [8,11,12]. Ниже ограничимся одним видом суррогатов, сгенерированных из временных рядов тех же систем, что использовались в качестве объектов моделирования, но без связи, и оценим статистическую значимость полученных значений показателя Р1 каждый раз с новыми начальными условиями и добавлением новой реализации динамического шума. Переходный процесс отсекаем.
2. Тестирование на эталонных примерах
В качестве аппроксимирующих функций в совместной и индивидуальной модели использовались полиномы произвольной степени со всевозможными перекрестными слагаемыми, как это делалось в [7]. Степень полинома варьировала от 1 до 6, размерность индивидуальной модели изменялась от 1 до 4, а размерность добавки -от 1 до 2. Такой подход нельзя назвать наилучшим во всех случаях, так как вполне возможно, что для тех или иных объектов предпочтительнее использовать иной базис, например, радиальные базисные функции [9]. Тем не менее, полиномиальный базис является самым распространенным и простым в использовании, требует меньшего объема данных, чем радиальные базисные или вейвлетные функции вследствие отсутствия локализации в фазовом пространстве. Рассматривалась ситуация, когда системы были связаны без запаздывания, то есть в формуле (2) т = 0.
Для оценки значимости использовались ансамбли из 100 рядов суррогатных данных и оценивалась 95% квантиль.
2.1. Связанные отображения окружности. В качестве объекта использовались два однонаправленно связанных отображения окружности (4) со связью вида k sin yn, поскольку переменная в отображении окружности имеет смысл фазы колебаний. Коэффициент связи k варьировал в пределах от 0.01 до 0.08 с шагом 0.002. В каждое отображение добавлялся нормально распределённый динамический шум с дисперсией, равной 0.001. Численным решением получались временные ряды длиной в 4096 точек, переходный процесс предварительно отсекался. Параметры для обоих отображений выбирались таким образом, чтобы при нулевой связи и отсутствии шума они соответствовали хаотическому поведению: Áx = 0.001, Áy = 0, Kx = 3.98, Ky = 3.88.
xn+1 — xn + Áx + Kx sin xn + k sin yn + Уп+1 = yn + Áy + Ky sin yn + Пп.
(4)
В состав отображения окружности входит непрерывная неполиномиальная функция 8т(жга). Поскольку данная функция может быть разложена в бесконечный степенной ряд, можно ожидать, что результат работы метода нелинейной причинности по Грейнджеру будет зависеть от степени полиномиальной функции, используемой в модели.
На рис. 1 представлены результаты применения метода причинности по Грейн-джеру - значения показателя улучшения прогноза в зависимости от уровня связи к. В качестве аппроксимирующих функций использовались полиномы различной степени от 1 до 6, размерность совместной модели равнялась 2 (03 = 1, Оа = 1).
Видно, что при использовании в качестве аппроксимирующей функции полиномов первой и второй степени серая и черная линии на графике зависимости Р1 (к) (определение наличия связи в верном и неверном направлении) принимают приблизительно одинаковые значения и лежат выше штриховой линии (95% суррогатного уровня).
Рис. 1. Зависимость показателя улучшения прогноза Р1 от величины коэффициента связи к при исследовании чувствительности метода на связанных отображениях окружности: непрерывной черной линией показаны значения Р1 при поиске связанности в верном направлении; серой линией - значения, полученные при поиске связи в заведомо ложном направлении (оценка воздействия системы X на У); штриховой линией - 95% суррогатный уровень
Рис. 2. Качественная проверка построенной модели: непрерывной черной линией показаны значения нелинейной функции отображения окружности, а штриховой серой - значения, полученные при её аппроксимации совместной моделью. Функция совместной модели д(х„,у„) рассчитывалась для всех пар значений (х„,у„) из экспериментального ряда, получался набор трёхмерных векторов (точек в трёхмерном пространстве) (х„,у„, д(хп, уп)), затем строилась проекция этого набора на плоскость (хп, д(х„,у„)), при этом точки сортировались по возрастанию хп
Рис. 3. Качественная проверка построенной модели. а - черной линией показаны значения наблюдаемого временного ряда отображения окружности (4); штриховой - значения временного ряда, восстановленного с помощью совместной модели вида (2); модель итерируется с начальными условиями, взятыми из исходного ряда (при построении предсказательной модели использовался полином степени 3 размерности 3). б - сплошной линией показаны значения наблюдаемого временного ряда воздействующей системы (ряд значений У)
Это говорит о том, что метод не может различить связь в верном и неверном направлении. Использование полиномов более высоких степеней позволяет диагностировать направленную связь в верном направлении, так как значения Р1 возрастают с ростом коэффициента связи к, оказываются значимыми и, более того, стремятся к 1. Это можно объяснить тем, что с увеличением степени полиномиальной функции происходит все более точная аппроксимация синусоидальных нелинейностей, включая нелинейность в связи.
В случаях, когда метод демонстрирует хорошую работоспособность, нелинейная функция исходной эталонной системы хорошо аппроксимируется входящим в предсказательную модель полиномом, что видно на рис. 2. Также можно сравнить динамику построенной совместной модели и исходной системы, временной ряд которой получен при значении коэффициента связи к = 0.05 (рис. 3, а). Для этого исходная временная реализация {хп} была наложена на реализацию, полученную путем итерирования реконструированного отображения (2) с начальными условиями, взятыми также из исходного ряда при использовании реализации {уп} (рис. 3, б) в качестве внешнего воздействия. Видно, что полученная в результате работы метода совместная модель достаточно хорошо качественно описывает динамику поведения исходной системы.
2.2. Связанные отображения Икеды. В качестве второго примера были взяты два отображения Икеды
Хп+1 = Ах + Бх[хп ео8(хП+иП+фх) -и,п вт(хП+иП+фх)] + к сов(у2п-хП) + 1п, Пп+1 = Бх [хп 8т(хп + ип + фх) + Ип е08(хп + ип + фх)], Уп+1 = АУ + БУ [Уп е°8(уп + ^п + фУ) - Уп й1п(уп + ^п + фУ)] + Пп, ^У [уп I
Уп+1 = БУ [уп 8т(уп + + фу) + Уп е°й(уп + ^п + фу
(5)
связанные однонаправленной нелинейной связью вида к еов(уп — хп), где у - первая переменная воздействующего отображения. Связь вносилась таким образом, потому что переменная в отображении Икеды имеет смысл фазы колебаний. Коэффициент связи к варьировал в пределах от 0.005 до 0.5 с шагом 0.005, к обоим отображениям добавлялся динамический шум и пп с дисперсией, равной 0.001.
Данная система сложнее для аппроксимации моделями вида (1), (2), поскольку кроме непрерывной неполиномиальной функции, разложимой в бесконечный степенной ряд, имеется также скрытая* переменная ип (ряд ип считался ненаблюдаемым, то есть не использовался при построении прогностической модели). Это делает невозможным точную или асимптотически точную аппроксимацию исходных эталонных уравнений, таким образом затрудняя решение задачи о поиске связанности.
От каждого отображения получали временной ряд длиною в 4096 точек, переходный процесс предварительно отсекался. По карте режимов выбирались параметры для обоих отображений такими, чтобы при нулевой связи и в отсутствие шума они соответствовали хаотическому поведению: Ах = 4, Вх = 0.2, фх = 0.5, Ау = 6, Ву = 0.13, фу = 0.4.
На графиках, представленных на рис. 4, показаны зависимости показателя улучшения прогноза от силы связи в случае использования в качестве аппроксимирующей модели полиномов разной степени (степень варьировала от 1 до 6), размерность совместной модели бралась равной 3 (Д5 = 2, Да = 1). Из графиков видно, что использование в качестве аппроксимирующей функции полинома первой степени не позволяет выявить направленную связь, так как значения показателя улучшения прогноза как в верном, так и в заведомо ложном направлении оказываются незначимыми, как и для случая связанных отображений окружности. Увеличение степени полинома позволяет выявить наличие связи. Значения показателя улучшения прогноза при поиске связи в ложном направлении не возрастают с ростом связи во всех случаях. Однако при больших степенях полинома, используемого в качестве аппроксимирующей функции, они оказываются значимыми, что может говорить о ненадёжности подхода, использующего слишком большие модели. Значения Р1 при поиске связи в верном направлении растут с ростом коэффициента связи к.
Отдельно следует отметить, что при использовании полиномов 5 и 6 степени показатели улучшения прогноза при малых значениях коэффициента связи принимают отрицательные значения. Теоретически это невозможно, но происходит на практике вследствие недостаточной точности вычислений. На наш взгляд, одним из факторов, приводящих к этому, является то, что использование полиномов большой степени при аппроксимации синусоидальной нелинейности приводит к большим ошибкам в определении минимума целевой функции. Используемая нами модель допускает применение линейного метода наименьших квадратов, то есть проблема не может быть обусловлена попаданием в локальный минимум и лежит, главным образом, в точности вычислений. В работах [15,16] показано, что при реализации метода наименьших квадратов с помощью рЯ-разложения, ортогонализации Грамма-Шмидта или в форме вращения Хаусхолдера ошибка в решении пропорциональна норме матрицы базисных функций. При использовании полиномов высоких степеней эта норма больше. Этот недостаток начинает сказываться, когда число коэффициентов оказывается слишком велико (сам размер матрицы возрастает), поэтому для совместной модели ошибка больше, чем для индивидуальной.
Поскольку не все переменные модели наблюдались, восстановить нелинейную функцию, как это было сделано для отображения окружности (см. рис. 2), невозможно. Потому были построены только временные ряды исходной системы при значении коэффициента связи к = 0.5 и реконструированной (рис. 5). Видно, что динами-
* Скрытыми называют переменные, входящие в уравнения модели, но недоступные измерению или наблюдению [13,14].
ка эталонной системы и предсказательной модели различны: неплохо аппроксимируя непосредственно следующее значение, модель в перспективе не повторяет качественно динамику объекта, а вместо хаотического ряда демонстрирует устойчивое положение равновесия, возмущаемое сигналом второго отображения. Тем не менее, метод нелинейной грейнджеровской причинности демонстрирует хорошие чувствительность и специфичность.
Рис. 4. Зависимость показателя улучшения прогноза Р1 от величины коэффициента связи к при исследовании чувствительности метода на связанных отображениях Икеды. Черной линией показаны значения Р1 при поиске связанности в верном направлении; серой линией - значения, полученные при поиске связи в заведомо ложном направлении (оценка воздействия системы X на У); штриховой линией - 95% суррогатный уровень
О 5 10 15 20 25 п
Рис. 5. Качественная проверка построенной модели: черной линией показаны значения наблюдаемого временного ряд отображения Ике-ды (5); штриховой серой - значения временного ряда, восстановленного с помощью совместной модели вида (2). Модель итерируется с начальными условиями, взятыми из исходного ряда. При построении предсказательной модели использовался полином степени 3 размерности 3
В приведённом на рис. 4 примере рассмотрен случай Ds = 2. Поскольку вектор состояния в данной ситуации нельзя восстановить точно, можно предполагать, что результаты могут различаться при различных Ds: например, ожидать, что оптимальные результаты могут быть достигнуты при размерности Ds = 5, согласно теореме Такенса [17]. Однако проведённые численные эксперименты показывают: размерности Ds = 2, как правило, достаточно для выявления существующей связи и признания связи в заведомо ложном направлении незначимой; при увеличении размерности результаты качественно не меняются.
2.3. Связанные отображения Заславского. Следующим объектом для исследования были выбраны два отображения Заславского [18] связанные односторонней линейной связью вида kyn, где y - первая переменная воздействующего отображения.
Xn+1 = (xn + Ax + Kx sin Xn + axu,n + kyn + bn) mod 2п, Un+i = aun + Kx sin Xn.
(6)
( yn+1 = (yn + Ay + Ky Sin yn + ayVn + Пп) mod 2п, l Vn+i = avn + Ky sin yn.
Коэффициент связи варьировал в пределах от 0.01 до 1.0 с шагом 0.01, к обоим отображениям добавлялся динамический шум с дисперсией 0.001. Нелинейные функции отображения Заславского являются более трудными для аппроксимации, чем функции отображений окружности и Икеды (наличие скрытой переменной и неполиномиальных функций, одна из которых является разрывной за счёт взятия остатка от деления на 2п). От каждого отображения получался временной ряд первой координаты длиною в 4096 точек, переходный процесс предварительно отсекался. Для того, чтобы получить хаотические режимы для систем X и Y, выбирались следующие значения параметров: Ax = п, Ay = 4, Kx = 4.59, Ky = 4.9, ax = 0.78, ay = 0.58.
Результаты работы метода причинности по Грейнджеру в случае использования в качестве аппроксимирующих функций полиномов 1...6 степени представлены на рис 6 (размерность совместной модели равнялась трем (Ds = 2, Da = 1)). В случае использования полинома 1 степени результаты оказываются значимы, в то время как для отображения Икеды подобная модель была неспособна обнаружить связь. На наш взгляд, это есть следствие введения линейной связи. Тем не менее, использование в качестве аппроксимирующих функций полиномов более высокой степени позволяет достигнуть лучшей специфичности - число неверных обнаружений связи в заведомо ложном направлении для них заметно меньше. Также при использовании полиномов степени 2 и выше показатель улучшения прогноза растет с увеличением значения коэффициента связи, в то время как при использовании полиномов степени 1 наблюдается уменьшение значения PI при больших коэффициентах k. Таким образом, даже при рассмотрении линейно связанных систем использование нелинейных моделей может давать некоторые преимущества.
Сравнивая временной ряд объекта, построенный при значении коэффициента связи k = 0.5, и ряд, сгенерированный восстановленной совместной моделью, можно видеть (рис. 7), что динамика модели и объекта существенно различна. Как и для
отображения Икеды, воспроизвести наблюдаемый режим качественно не удается. Это не удивительно, поскольку вектор состояния наблюдался не полностью, а тригонометрические нелинейные функции с разрывами аппроксимировались полиномами низких степеней. Несмотря на это метод причинности по Грейнджеру демонстрирует разумные результаты.
Рис. 6. Зависимость показателя улучшения прогноза Р1 от величины коэффициента связи к при исследовании чувствительности метода на связанных отображениях Заславского. Сплошной черной линией показаны значения Р1 при поиске связанности в верном направлении; серой линией - значения, полученные при поиске связи в заведомо ложном направлении (оценка воздействия системы X на У); штриховой линией - 95% суррогатный уровень
Рис. 7. Качественная проверка построенной модели: черной линией показаны значения наблюдаемого временного ряда отображения Заславского (6); штриховой серой - значения временного ряда, восстановленного с помощью совместной модели вида (2). Модель итерируется с начальными условиями, взятыми из исходного ряда. При построении предсказательной модели использовался полином степени 3 размерности 3
2.4. Связанные системы Лоренца. В предыдущих разделах метод причинности по Грейнджеру применяли для определения связанности по рядам эталонных систем с дискретным временем. Рассмотрим в качестве объекта две связанные однонаправленной связью системы с непрерывным временем - системы Лоренца
Ж1 = Oi(x2 - Xi), Ж2 = Xi(ri - Жэ) - X2, жэ = Ж1Ж2 - biX3 + k/э,
(7)
y/1 = o2(/2 - yi),
y/2 = /i(r2 - Уэ) - /2, //э = yi/2 - b2/3-
Следует отметить, что в принципе такая постановка задачи не слишком отличается от ранее рассмотренных, поскольку аппроксимация дифференциальных уравнений разностной схемой возможна (требуется проверить также условие сходимости). Однако в рассмотренном случае две из трёх переменных каждой подсистемы в (7) считались ненаблюдаемыми, что существенно усложняет построение модели. Ситуация аналогична рассмотренной для двух связанных отображений Икеды. Более того, дополнительно присутствует специфика, связанная с тем, что временные ряды получены от системы с непрерывным временем: большие времена корреляции, спектр с выраженными максимумами, большое расстояние между экстремумами в реализации. Эти особенности отличают данную ситуацию от ранее рассмотренных случаев.
Линейная аддитивная связь вида k/э вносилась в уравнение для третьей координаты Жэ, коэффициент связи k варьировал в пределах от 0.01 до 1 с шагом 0.01, шаг интегрирования 0.001. От каждой системы получали временные ряды длиною в 20480 точек на аттракторе. Хаотические режимы реализовывались в системах с параметрами: oi = 10, о2 = 10.5, ri = 30, r2 = 31, bi = 8/3, b2 = 8/3.
Для графиков, представленных на рис. 8, в качестве моделей использовались полиномы степени 1 и 2 с индивидуальной и совместной размерностями, равными 1 (а, б) и 2 (в, г). Лаг l = 1.0. В случае использования в качестве линейной модели полинома степени 1 (рис. 8, а, б) результаты оказываются значимы только при больших размерностях (см. рис. 8, б) и при значениях коэффициента связи k ~ 0.8. Нелинейная модель (см. рис. 8, в, г) с полиномом степени 2 в качестве аппроксимирующей функции позволяет выявить связь (см. рис. 8, г), начиная с вдвое меньших значений коэффициента связи k ~ 0.5, а значения показателя улучшения прогноза PI оказываются больше, чем при использовании линейной модели.
При поиске связи в заведомо ложном направлении присутствует небольшой рост значений PI, однако результаты оказываются незначимы, с точки зрения суррогатов, для всех рассмотренных коэффициентов связи.
На рис. 9 представлено сравнение исходного временного ряда, построенного при значении коэффициента связи k = 0.6, и ряда, полученного с помощью совместной модели. Видно, что уравнения (1), (2) полностью не описывают динамику системы. Однако метод нелинейной причинности по Грейнджеру и в этом случае позволяет получить значимые результаты в заведомо верном направлении и незначимые - в заведомо ложном, то есть проходит проверку как на чувствительность, так и на специфичность, начиная с некоторого значения коэффициента связи.
Рис. 8. Зависимость показателя улучшения прогноза Р1 от величины коэффициента связи к при исследовании чувствительности метода на связанных однонаправленной связью системах Лоренца при использовании: а, б - линейной модели (размерность совместной модели 2 и 4, соответственно); в, г -нелинейной модели со степенью полинома 2 (размерность совместной модели 2 и 4, соответственно). Непрерывной черной линией показаны значения Р1 при поиске связанности в верном направлении; серой линией - значения, полученные при поиске связи в заведомо ложном направлении (оценка воздействия системы X на У); штриховой линией - 95% суррогатный уровень
10.0-1-,-,-,-,-,-,-,-,-,
0 1000 2000 3000 4000 п
Рис. 9. Качественная проверка сходства временных рядов: наблюдаемый временной ряд третьей координаты системы Лоренца (черная линия) и временной ряд, восстановленный совместной моделью вида (серая штриховая линия); модель интегрируется с начальными условиями, взятыми из исходного ряда. При построении предсказательной модели использовался полином степени 2, размерности 4
Заключение
На примере ряда последовательно усложняющихся эталонных однонаправлен-но связанных систем (системы, содержащие неполиномиальные функции, скрытые переменные, разрывные функции, описываемые дифференциальными уравнениями) изучалась работоспособность метода причинности по Грейнджеру. В качестве аппроксимирующих функций использовались полиномы общего вида. Метод в таком виде уже неоднократно применялся для различных реальных объектов, в частности, нейрофизиологических и климатологических. Целью исследования было уточнение границ его применения: насколько несоответствие структуры прогностической модели структуре породившего временные ряды объекта существенно сказывается на чувствительности и специфичности метода.
На ряде численных примеров было показано следующее.
• Работоспособность метода зависит от качества аппроксимации как собственных нелинейных функций модели, так и внешнего воздействия. Слишком сильное упрощение этих функций в предсказательной модели может привести к падению специфичности: метод перестаёт различать связь в верном и ложном направлениях. Лучшая аппроксимация нелинейных функций приводит к возможности выявления связи при меньших ее величинах.
• При достаточно высоких степенях и размерностях полинома метод верно диагностирует связь даже в случаях, когда прогностическая модель существенно отличается от эталонных уравнений из-за наличия скрытых переменных. Работоспособность метода сохраняется для временных рядов с различными статистическими характеристиками, типичными как для систем с дискретным временем, так и с непрерывным. Анализ аттрактора реконструированной прогностической модели показывает, что ее временной ряд имеет поведение принципиально непохожее на поведение ряда эталонной системы: вместо хаотического режима демонстрирует устойчивое равновесие или периодический режим, возмущаемый воздействием другой подсистемы. Таким образом, оказывается, что работоспособность метода напрямую не зависит от способности прогностической модели воспроизводить режим поведения наблюдаемой системы, хотя и может быть с нею связана; детальное прояснение этого момента требует отдельного исследования.
• Большое число коэффициентов полиномиальной модели приводит к нежелательным последствиям. Во-первых, предсказательная модель начинает описывать не процесс, а конкретные особенности измеренного ряда, в том числе, конкретную реализацию шума. Это является одной из причин того, что метод начинает показывать значимую связь в заведомо ложном направлении. Во-вторых, возрастают ошибки определения коэффициентов, поэтому минимум целевой функции определяется неточно и, как следствие, рассчитанное улучшение прогноза может принимать произвольные, в том числе, отрицательные значения.
По результатам работы можно сделать общие выводы и рекомендации.
• Использование нелинейных моделей в целом повышают специфичность и чувствительность метода.
• Можно рассчитывать на работоспособность подхода даже в случае, когда добиться качественного воспроизведения динамики наблюдаемого процесса не удаётся, то есть добиваться соответствия топологических характеристик аттракторов объекта и модели не всегда необходимо, хотя данный вопрос заслуживает отдельного исследования.
• Нужно следить за точностью вычислений и числом коэффициентов полиномиальной модели, поскольку возможность переобучения модели и рост влияния численных ошибок могут привести к неверным и даже нефизичным выводам. Одним из индикаторов такого поведения является получение отрицательных значений улучшения прогноза Р1, теоретически невозможных.
Авторы выражают благодарность профессору Б.П. Безручко, с которым неоднократно обсуждались результаты работы и чьи замечания и предложения помогли существенно улучшить ее.
Работа выполнена при поддержке грантов РФФИ: № 11-02-00599 и № 12-02-00377.
Библиографический список
1. Granger C.W.J. Investigating causal relations by econometric models and cross-spectral methods // Econometrica. 1969. Vol. 37, № 3. P. 424.
2. Andrea Brovelli, Mingzhou Ding, Anders Ledberg, Yonghong Chen, Richard Naka-mura, and Steven L. Bressler. Beta oscillations in a large-scale sensorimotor cortical network: Directional influences revealed by Granger causality // PNAS. 2004. Vol. 101. P. 9849.
3. L.A. Baccala, K. Sameshima, G. Ballester, A.C. Do Valle and C. Timo-Laria. Studing the interactions between brain structures via directed coherence and Granger causality // Applied sig. processing. 1998. Vol. 5. P. 40.
4. P. Tass, D. Smirnov, A. Karavaev, U. Barnikol, T. Barnikol, I. Adamchic, C. Hauptmann, N. Pawelcyzk, M. Maarouf, V. Sturm, H.-J. Freund, and B. Bezruchko. The causal relationship between subcortical local field potential oscillations and parkinsonian resting tremor // J. Neural Eng. 2010. Vol. 7. 016009.
5. И.И. Мохов, Д.А. Смирнов, П.И. Наконечный, С.С. Козленко, Ю. Куртс. Оценка взаимного воздействия Эль-Ниньо - Южного колебания и Индийского муссона // в «Современные проблемы динамики океана и атмосферы» / Ред. А.В. Фролов и Ю.Д. Реснянский. М.: ТРИАДА ЛТД, 2010. С. 251.
6. С.С. Козленко, И.И. Мохов, Д.А. Смирнов. Анализ причинно-следственных связей между Эль-Ниньо в Тихом океане и его аналогом в экваториальной Атлантике // Известия РАН. Физика атмосферы и океана. 2009. Т. 42, № 6. C. 754.
7. Yonghong Chen, Govindan Rangarajan, Jianfeng Feng, Mingzhou Ding. Analyzing Multiple Nonlinear Time Series with Extended Granger Causality // Physics Letters A. Vol. 324, Issue 1. P. 26.
8. И.В. Сысоев, А.С. Караваев, П.И. Наконечный. Роль нелинейности модели в диагностике связей при патологическом треморе методом грейнджеровской причинности // Изв. вузов. Прикладная нелинейная динамика. 2010. Т. 18, № 4. С. 81.
9. Marinazzo Daniele, Pellicoro Mario, and Stramaglia Sebastiano. Nonlinear parametric model for Granger causality of time series// Phys. Rev. E. 2006. Vol. 73. 066216.
10. Смирнов Д.А. Выявление нелинейных связей между стохастическими осцилляторами по временным рядам // Известия вузов. Прикладная нелинейная динамика, 2010. Т. 18, в. 2. С. 16.
11. Schreiber T. and Schmitz A. Surrogate time series // Physica D. 2000. Vol. 142. 346.
12. Kevin T. Dolan and Alexander Neiman. Surrogate analysis of multichannel data with frequency dependant time lag // Physical review E. Vol 65. 026108.
13. Baake E., Baake M., Bock H.G., and Briggs K.M. Fitting ordinary differential equations to chaotic data // Phys. Rev. A. 1992. Vol. 45, № 8. P. 5524.
14. Boris P. Bezruchko, Dmitry A. Smirnov and Ilya V. Sysoev. Identification of chaotic systems with hidden variables (modified Bock's algorithm) // Chaos, Solitons & Fractals. 2006. Vol. 29. P. 82.
15. Bjork A. Solving Linear Squares Problem by Gram-Schmidt Orthogona lization // Math. Copm. 1976. Vol. 20. P. 325.
16. Дж. Голуб, Ч. Ван Лоун. Матричные вычисления: Пер. с англ. М.: Мир, 1999. 548 с.
17. Takens F. Detecting strange attractors in turbulence // Lecture Notes in Math. 1981. Vol. 898. P. 366.
18. Заславский Г.М., Сагдеев Р.З., Усиков Д.А., Черников А.А. Слабый хаос и квазирегулярные структуры. М.: Наука, 1991. 236 с.
Саратовский госуниверситет Поступила в редакцию 7.11.2011
СФ ИРЭ им. В.А. Котельникова РАН После доработки 28.03.2013
INFLUENCE OF THE CHOICE OF THE MODEL STRUCTURE FOR WORKING CAPACITY OF NONLINEAR GRANGER CAUSALITY APPROACH
Maxim V. Kornilov, Ilya V. Sysoev
Currently, the method of nonlinear Granger causality is actively used in many applications in medicine, biology, physics, to identify the coupling between objects from the records of their oscillations (time series) using forecasting models. In this paper the impact of choosing the model structure on the method performance is investigated. The possibility of obtaining reliable estimates of coupling is numerically demonstrated, even if the structure of the constructed forecasting model differs from that of the reference system.
Keywords: The method of nonlinear Granger causality, the reconstruction of the time series, nonlinear dynamical systems.
Корнилов Максим Вячеславович - родился в 1988 году в Саратове, окончил Лицей математики и информатики (2005), механико-математический факультет Саратовского государственного университета имени Н.Г. Чернышевского (2009, присуждена степень бакалавра математики), факультет нано- и биомедицинских технологий (2011, присуждена степень магистра техники и технологии по направлению биомедицинская инженерия). В настоящее время аспирант кафедры динамического моделирования и биомедицинской инженерии. Научные интересы: анализ временных рядов, математическое моделирование биологических процессов, автоматическое распознавание изображений, компьютерное зрение.
410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: [email protected]
Сысоев Илья Вячеславович - родился в Саратове (1983). Окончил Лицей прикладных наук (1999) и факультет нелинейных процессов СГУ (2004). Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2007). Работал на кафедре электроники, колебаний и волн (2005-2007). В настоящее время - доцент базовой кафедры динамического моделирования и биомедицинской инженерии. Научные интересы - исследование сигналов биологической природы методами нелинейной динамики, исследование эффективности и модернизация подходов к анализу сигналов. Автор более 40 публикаций.
410019 Саратов, ул. Зелёная, 38
Саратовский филиал Института радиотехники и электроники им. В.А. Котельникова РАН E-mail: [email protected]
