Свойства одного алгоритма оценки параметров нестационарных пуассоновских процессов Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 519.22

СВОЙСТВА ОДНОГО АЛГОРИТМА ОЦЕНКИ ПАРАМЕТРОВ НЕСТАЦИОНАРНЫХ ПУАССОНОВСКИХ ПРОЦЕССОВ *

Г.М. Водинчар1,2

1 Камчатский государственный университет имени Витуса Беринга, 683032, Петропавловск-

Камчатский, ул. Пограничная, 4

2 Институт космофизических исследований и распространения радиоволн ДВО РАН,684034,

Камчатский край, с. Паратунка, ул. Мирная, 7

E-mail: [email protected]

Рассмотрены свойства алгоритма оценки параметров гармонических составляющих интенсивности нестационарного пуассоновского процесса. Установлена хорошая вычислительная устойчивость алгоритма, оценена скорость сходимости оценок. Рассмотрена возможность применения алгоритма к анализу приливного отклика в импульсном электромагнитном излучении Земли.

Ключевые слова: метод наименьших квадратов, нестационарный пуассонов-ский процесс, импульсное электромагнитное излучение Земли.

(с) Водинчар Г.М., 2012 MATHEMATICAL SIMULATION

MSC 62M10

CHARACTERISTICS OF AN ALGORITM FOR NON-STATIONARY POISSON PROCESS ESTIMATION

G.M. Vodinchar1, 2

1 Vitus Bering Kamchatka State University, 683032, Petropavlovsk-Kamchatsky, Pogranichnaya

st., 4, Russia

2 Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch

Russian Academy of Sciences, 684034, Kamchatsky Kray, Paratunka, Mirnaya st., 7,

Russia

E-mail: [email protected]

The properties of the algorithm estimates of parameter of harmonic intensity for nonstationary Poisson process. The computational stability of the algorithm is proved.

Key words: least squares method, non-stationary Poisson process, electromagnetic pulsed terrestrial emission.

(c) Vodinchar G.M., 2012

*Работа выполнена при поддержке Минобрнауки России по "Программе стратегического развития КамГУ им. Витуса Беринга"на 2012-2016 гг.

Введение

В настоящей работе исследуются некоторые вычислительные аспекты разработанного автором в [1] алгоритма анализа временных рядов, описывающего один класс нестационарных пуассоновских процессов. Речь идет о процессах, интенсивность которых меняется во временем по полигармоническому закону, с априори известным частотным составом. Такие модели возникают при описании точечных процессов в системах, подверженных внешним воздействиям, спектр которых известен. В геофизических системах подобным воздействием являются волны приливного потенциала, состав которого известен детально [2]. Если величина отклика геофизического поля на воздействие в среднем близка к фоновым шумам, то усиление этого отклика является индикатором приближения геофизической среды среды к критическому состоянию [3]. Алгоритмы обнаружения приливного отклика в импульсном электромагнитном излучении Земли, их статистические свойства, результаты тестирования и результаты применения к реальным данным описаны в [1]. Одним из этапов алгоритмов является расчет точечных и интервальных оценок для амплитуд и фаз гармонических составляющих в потоке импульсов. В данной статье рассмотрены вопросы вычислительной обусловленности и скорости сходимости оценок.

Постановка задачи оценивания и описание алгоритма

Рассматривается нестационарный пуассоновский процесс ^ (?), интенсивность которого является суммой постоянной составляющей и нескольких гармоник. Реализация процесса описывается временным рядом х - числом скачков процесса в интервале времени (?;? + 1]. Тогда компоненты ряда - независимые пуассоновские случайные величины с интенсивностью

Частоты предполагаются известными, амплитуды Ак и Вк подлежат оценке во временном окне из 2N +1 отсчетов. Использование нечетного числа отсчетов обусловлено вычислительными удобствами. Всегда можно также считать, что — N < ? < N.

Формально представляя члены ряда хг в виде хг = , где £г = хг — , полу-

чаем задачу оценки параметров тренда ^ на фоне нестационарного белого шума £{. Шум, очевидно, негауссовский, имеет нулевое среднее и дисперсию ^. Таким образом, распределение шума зависит от оцениваемых параметров. Для решения задачи применяется метод наименьших квадратов (МНК), т.е. минимизация функции

N

^ (хг — Х1 )2 по амплитудам Ак, Вк, Ао. Негауссововсть и нестационарность шума,

зависимость его дисперсии от неизвестных параметров не позволяют использовать для характеризации оценок общую теорию МНК-оценивания [5].

Для сокращения объема вычислений при обработке ряда в скользящем временном окне удобно перейти к базисным функциям с нулевым средним по времени. Вводим обозначения

Xt = Ao + £ (AkCosOkt + BkSino>kt)

k=1

(1)

N

akt = cGs^kt

(2N + 1)sin0.5ak7

sin(N + 0.5) a>k

ßkt = sino>kt, C = A0 +

Тогда выражение (1) примет вид

^г = С + £ <Акакг + Вкркг) (2)

к=1

Выражая синусы и косинусы через комплексные экспоненты по формулам Эйлера легко установить непосредственным вычислением, что базисные функции акх

N

и рр ортогональны относительно скалярного произведения {ак, в]) = £ акгрр при

г=—N

любых к и ] и имеют нулевое среднее по времени для — N < г < N. В работе [1] показано, что МНК-оценки А*, В* и С* параметров Ак, Вк и С соответственно являются несмещенными, сходятся в среднеквадратическом смысле к истинным значениям и асимптотически нормальны.

Вычислительная обусловленность МНК-оценивания

Ввиду вышеуказанной частичной ортогональности базисных функций МНК-система для оценок амплитуд распадается на три независимые подсистемы

N р р

^ + 1)С = £ Хг, £ (а], а^А] = (х, ар), £ (рр, р1)Вр = (х, ар), I = 1,..., р. (3)

t=—N ]=1 ]=1

Из этих формул видно, что с вычислительной устойчивостью С* трудностей нет, а для А* и В* возникает вопрос об обусловленности матриц \\(ар, аг-)\| и \\(рр,рг-)||. Непосредственным вычислением можно получить следующие формулы для элементов этих матриц:

(а,а) = N + 1/2 +sln<2N +1)Ю sin2(N+°,5)Ю ■

281пюг- <2N + 1^т2°,

sin<N + 0,5)<Юг- — Юр) sin<N + °, 5)<Юг- + Юр) sin<N + °, 5)Юг^т^ + 0.5) Юр

а, ар 2sin°, 5<Юг- — Юр) + 2sin°, 5<Юг- + Юр) <2N + 1)sin°, 5Юг^т°, 5юр ’

„,.р.) =N+,/2— («

р sin<N + 0 ,5)<юг- — Юр) sin<N + °.5)<юг- + Юр) , = ,

, р 2sin°, 5<Юг- — Юр) 2sin°, 5<Юг- + Юр) , р

Из этих выражений видно, что для достаточно длинного временного окна в матрицах \\(ар, аг-)|| и \(рр,рг-)\| будет диагональное преобладание, что гарантирует их хорошую обусловленность и устойчивое вычисление МНК-оценок. Ясно также, что величина окна, при котором гарантировано диагональное преобладание определяется близостью различных частот, поскольку при фиксированном N внедиагональные <1,^-элементы при неограниченном сближении частот Юг- и Юр растут как N, т.е. сравнимы с диагональными элементами.

При разделении гармоник при обработке временных рядов величина окна должна превосходить период биения гармоник самых близких частот. Рассмотрим более этот вопрос на примере ряда интенсивности импульсного электромагнитного излучения

Таблица 1

Периоды основных приливных волн и периоды биений

k Волна Гармоника Период Тк, мин Частота Юк, рад/мин Период биений для Юк и Юк— 1, сут

1 Ol 1 1549,160455 0,004055865 -

2 Kl 1 1436,068141 0,00437527 13,66079164

3 Ol 2 774,5802275 0,00811173 1,167769287

4 M2 1 745,236078 0,008431134 13,66080398

5 S2 1 720 0,008726646 14,76529114

6 Kl 2 718,0340705 0,008750539 182,6194862

7 Ol 3 516,3868183 0,012167594 1,276925027

В Kl 3 478,6893803 0,013125809 4,553597214

9 Ol 4 387,2901138 0,016223459 1,408591297

10 M2 2 372,618039 0,016862268 6,830401989

11 S2 2 360 0,017453293 7,382645572

12 Kl 4 359,0170353 0,017501078 91,30974311

13 M2 3 248,412026 0,025293402 0,559951464

14 S2 3 240 0,026179939 4,921763715

15 M2 4 186,3090195 0.033724536 0.578337422

16 S2 4 180 0,034906585 3,691322786

Земли, для анализа которого изначально разрабатывался алгоритм [1]. В этом ряде данных шаг по времени составлял 1 мин. Исследовался отклик на наиболее сильные приливные волны 01, Кь М2, $2 [2] и их высшие гармоники до 4 включительно. Периоды этих волн и периоды биения наиболее близких по частотам гармоник приведены в Таблице 1. Видно, что периоды биения составляют от полусуток до полугода.

Непосредственные вычисления показывают, что для всех этих случаев, если длина окна превосходит период биения, то матрицы МНК-систем имеют диагональное преобладание. Поэтому при вычислении оценок можно использовать обычные конечные методы линейной алгебры, без использования регуляризующих процедур.

Скорость сходимости МНК-оценок

Обозначим матрицы, обратные к ||(а/,аг-)\| и \(рр,рг-)\| через и и V, соответственно. В [1] показано, что элементы этих матриц бесконечно малы при N ^ <*>.

Скорость среднеквадратической сходимости МНК-оценок характеризуется их дисперсиями. Теоретические оценки этой скорости даются следующим утверждением.

Теорема 1. БС* является функцией порядка 1/N при N ^ <*>. Для ПА* и ПВ* такая функция является асимптотически верхней оценкой.

Доказательство. Прежде всего отметим, что функция Хг ограничена сверху некоторым положительным числом Л, не зависящим от N. Из первой из систем (3) видно,

N

N

N

что С* = 2мГГ £ X ’ откуда БС* =

2ІУ + 1 ¿=_ N

(2Ñ+ÍF Е, Dx'

(2N + 1)2 t=EN^

л

2N +1’

что и завершает доказательство для С*.

Получим теперь верхнюю оценку модулей диагональных элементов матрицы и. Поскольку матрица ||(оу, аг-)|| симметрична и положительно определена, она обладает ортонормированной системой собственных вектор-столбцов Як положительных собственных значений Ик. Тогда матрица и допускает разложение [6]

U = H^j, аi)

11 s1... sp

ß1 Vp

[s1... sp]T

С учетом того, что все компоненты векторов Як не превосходят по модулю единицу,

1,1 1 р р

из этого разложения видно, что |ий| <-1--1---<----- —т- < —, где то - нижняя

И Ир шт|Дк} то

граница области Гершгорина для матрицы ||(aj, аг-)||.

Напомним, что область Гершгорина (содержащая все собственные значения матрицы) в комплексной плоскости является объединением кругов, центры которых -диагональные элементы, а радиусы - суммы модулей недиагональных элементов соответствующих строк [7]. Из формул (4) видно, что недиагональные элементы матрицы Ц(а, аг-)|| ограничены, а диагональные являются суммой N и ограниченной величины. Тогда область Гершгорина при N ^ <^, оставаясь ограниченной, удаляется в бесконечность, точнее говоря существует такая положительная константа К, что то < N — К, откуда |иг-г-| <

p

N - К

Теперь выполним следующую цепочку преобразований

N

N

DA* = D E Ukj E аjtXt = D E xt E Ukj<

аjt =

j=1

N

-N j=1

N

EN

t= N

E Ukj аjt

Dxt =

N

N

E ( E Ukj ар\ ^ — л E ( E Ukj

аjtі =

t= N j=

j=1

p

t=-N \j=1 p

Np

л ^E Ukjukiаjtай

t=-N j,i=1

N p p

л ^E UkjUki аjtаit л ^E Ukjuki{аj, аі) л Uki&ki Лukk.

j,i=1 t=-N j,i=1 i=1

лp

N — K.

Доказательство теоремы

Из этих соотношений видно, что Ыкк > 0 и БА* <

для БВ* полностью аналогично. □

Полученная в теореме оценка дает верхнюю границу среднеквадратической ошибки в вычислении БА* и БВ*. Входящая в формулу оценки константа К для конкретного набора частот и длины окна может оказаться достаточно большой и сравнимой с N .В связи с этим интересно установить в вычислительных экспериментах насколько близки могут быть |Ыц| и |Ыц| к своей теоретической верхней границе ———, т.е.

N — К

насколько завышена эта теоретическая оценка.

Были проведены вычисления чисел Ы = тт{|ыг7|} и V = тІп{|Уй|} для первых гармоник приливных волн 0і, Кі, М2, $2 (общее число гармоник — = 4) при N от 104 до 5 ■ 104 с с шагом 103, что соответствует величине временного окна от 14 до 70 сут. Отметим, что минимальное значение N, необходимые для разделения гармоник вышеуказанных приливных волн составляет 10630.

1

2

t

t

2

2

Полученная зависимости и от N в двойном логарифмическом масштабе и прямая регрессии приведены на рис. 1. Аналогичная зависимость V от N отличается от нее настолько мало, что при совмещении на одном рисунке расчетные точки и прямые регрессии визуально совпадают.

Видно, что обе зависимости имеют четко степенной характер. Наклон аппроксимирующей прямой для и^) составляет —1,00007, а для и^) этот наклон равен — 1,00003. Таким образом, теоретическая оценка убывания погрешности по закону 1/N практически совпадает с расчетной.

1пи

\

ч

\

'• InJV

9 9,2 9,4 9,6 9,8 10 10,2 10,4 10,6 10,8 11

Зависимость и от N в двойном логарифмическом масштабе и прямая регрессии.

Заключение

В работе проведено исследование свойств алгоритма МНК-оценивания параметров гармонических компонент интенсивности нестационарного пуассоновского процесса, описываемого временным рядом числа скачков в единицу времени. Получены теоретические оценки обусловленности матрицы МНК-системы и среднеквадратической скорости сходимости оценок к точным значениям. Продемонстрирована возможность применения алгоритма к выявлению приливного отклика в импульсных геофизических сигналах на примере импульсного электромагнитного излучения Земли.

Библиографический список

1. Водинчар Г.М. Методика оценивания некоторых параметров импульсного электромагнитного излучения Земли // Вычислительные технологии. 2002. № 5. T. 7. C. 29-35.

2. Мельхиор П. Земные приливы. М.: Мир, 1968.

3. Любушин А.А. Анализ данных систем геофизического и экологического мониторинга. М.: Наука, 2007.

4. Водинчар Г.М. Алгоритмы и программы оценивания параметров гармонических составляющих временных рядов пуассоновского характера: дис. ... канд. физ.-мат. наук. Петропавловск-Камчатский: КГПУ, 2002.

5. Айвазян С.А., Енюков И.С., Мешалкин Л.Д. Прикладная статистика: исследование зависимостей. М.: Финансы и статистика, 1985.

6. Марпл-мл С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990.

7. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1967.

Поступила в редакцию / Original article submitted: 03.08.2012