CC BY
9
ISSN 1992-6502 (Print)_
2020. Т. 24, № 1 (87). С. 93-102
Вестник УГАТУ
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 519.676
Приближенное решение задачи оптимальной фильтрации по МАП-критерию с применением частичных сумм ряда Эджворта
к. а. рЫБАКОВ
rkoffice@mail.ru
ФГБОУ ВО «Московский авиационный институт (национальный исследовательский университет)» (МАИ)
Поступила в редакцию 20.01.2020
Аннотация. На основе непрерывного фильтра частиц и аппроксимации плотности вероятности частичными суммами ряда Эджворта сформирован статистический алгоритм приближенного решения задачи оптимальной фильтрации по критерию максимума апостериорной плотности вероятности. В предложенном алгоритме мода апостериорного распределения для каждой координаты ненаблюдаемого вектора состояния динамической системы неявно выражается через апостериорные центральные моменты до третьего, четвертого, пятого или шестого порядков включительно. Приводится модельный пример.
Ключевые слова: максимум апостериорной плотности вероятности; мода распределения; полиномы Эрмита; оптимальное оценивание; оптимальная фильтрация; ряд Эджворта; статистический алгоритм; стохастическое дифференциальное уравнение; стохастическая система; фильтр частиц.
ВВЕДЕНИЕ
Оптимальная фильтрация сигналов - одна из важнейших задач, возникающая при проектировании и внедрении систем управления динамическими системами, которые действуют в условиях случайных возмущений и помех [ 1, 2].
В статье предлагается статистический алгоритм фильтрации сигналов в стохастических системах наблюдения и оценивания с непрерывным временем. Предполагается, что модель объекта наблюдения и модель измерительной системы задаются стохастическими дифференциальными уравнениями, оптимальность оценки определяется по критерию максимума апостериорной плотности вероятности ненаблюдаемого вектора состояния (МАП-критерий [3]). В основе предлагаемого алгоритма лежит непрерывный фильтр частиц [4] и аппроксимация апостериорной плотности вероятности ча-
стичными суммами ряда Эджворта [5-7]. Такое представление используется для маргинальных плотностей, оно удобно тем, что выражается через центральные моменты апостериорного распределения или через соответствующие кумулянты [8], а их можно оценить, моделируя объект наблюдения и весовые коэффициенты, возникающие в алгоритме фильтра частиц [7].
Статья является продолжением работы [7], в которой частичная сумма ряда Эджворта определялась моментами до 3-го порядка включительно, что было обусловлено удобством применения аналитических формул решения алгебраического уравнения 4-й степени для приближенного нахождения моды апостериорного распределения. В этой статье учитываются моменты до 6-го порядка включительно. Отметим, что представление плотности вероятности рядом Эджворта или его частичной суммой при-
Работа поддержана грантом РФФИ 17-08-00530-а.
менялось как в задачах анализа стохастических систем, так и в задачах фильтрации [9-11]. Использование моментов до 4-го или 6-го порядков включительно предлагалось в [11, 12], но сами моменты требовалось находить как решение усеченной системы дифференциальных уравнений, полученной с помощью специальной процедуры замыкания и метода статистической аппроксимации нелинейностей.
В данном случае сформирован статистический алгоритм, в котором моменты оцениваются методом Монте-Карло. Такой подход оправдан, поскольку современные вычислительные средства позволяют моделировать сложные многомерные стохастические системы в реальном времени [13].
Отметим, что фильтр частиц предполагает моделирование ансамбля траекторий объекта наблюдения и соответствующих весовых коэффициентов, по которым можно оценить не только моменты распределения, но и апостериорную плотность вероятности, например с помощью гистограммы или ядерной оценки [14]. Это дает возможность решить задачу оптимальной фильтрации по различным критериям сколь угодно точно. Однако оценка апостериорной плотности вероятности по сравнению с оценкой моментов требует увеличения числа моделируемых траекторий в ансамбле. В первую очередь это связано с особенностями согласованного выбора параметров при оценивании (шага численного интегрирования, объема выборки и др.) [15]. В приложении к задаче оптимальной фильтрации согласованный выбор параметров обсуждался в работе [7]. Кроме того, это более трудоемкая процедура, реализовать которую в реальном времени затруднительно. Таким образом, приближенное решение задачи оптимальной фильтрации по МАП-критерию с применением частичных сумм ряда Эджворта может служить компромиссом между точностью оценивания и вычислительными затратами.
ПРИБЛИЖЕННОЕ НАХОЖДЕНИЕ МОДЫ РАСПРЕДЕЛЕНИЯ
Пусть X - скалярная случайная величина, для которой существуют моменты всех порядков. В частности, Ми Б - ее матема-
тическое ожидание и дисперсия соответственно: М = ЕХ, Б = Е(Х- М)2, Е - оператор математического ожидания. Определим центрированную и нормированную случайную величину
5=Х=М. а = 4Б.
а
Плотность вероятности ф(^) случайной величины можно представить рядом Эджворта [8, 16]: ф® = ф® х
1 Д3
1 ( Д4
1 + ^ ^ Н з(0 + - 3 I Н 4® +
3! а
(3)
4!I а4
(4)
10 ( д3
+ ыI Н6®+5,105- 1Н5)+
(4)
(5)
35 д3 ( д
280 ( д3
ЯН ^-31Н7(^)^1 н9+
7! а
9!
(5)
(5)
1
+ — 6!
(
-15 -10 (I + 30
ч а6 а4 и31 ,
(6)
ч2
Л
Н 6(^)+
56 ( д, -10 дЛ 35 (д^ - 3
8! а31а5 10а31+ 8! |а4 3
(6)
+^ & I & - 3) Н^+
(6)
Н Д) +
+^ &1НИ®+• ^
(6)
где ф(£) - плотность вероятности случайной величины, имеющей стандартное нормальное распределение (с нулевым математическим ожиданием и единичной дисперсией), Нк($ - полином Эрмита степени к:
Ф® =
1
л/2л
й
Ф к^ = Ф(^) = (-1)к Ф(^)Нк
и др - центральные моменты порядка р случайной величины X, т. е. др = Е(Х - М)р. При учете центральных моментов др до порядка р включительно в частичной сумме ряда
Эджворта необходимо оставлять слагаемые, отмеченные как (г) для всех г <р.
На основе этого разложения в [8] была предложена следующая аппроксимация моды распределения с плотностью вероятности ф(£) при учете центрального момента ц3 (асимметрия Шарлье):
1 М 2 а3
что соответствует аппроксимации моды распределения случайной величины X:
х * = м -1М-.
2 а2
(1)
В [5-7] дополнительно была использована другая аппроксимация, а именно точное значение с, , на котором частичная сумма ряда Эджворта при учете центрального момента ц3 достигает максимального значения. При таком подходе, во-первых, есть возможность провести сравнение с приближением (1), а во-вторых, учет центрального момента ц3 приводит к необходимости решения алгебраического уравнения 4-й степени и можно воспользоваться методами, позволяющими найти точно его корни. Учет центральных моментов более высокого порядка приводит к необходимости решения алгебраических уравнений, степень которых выше 4, т. е. требуется применять приближенные методы нахождения корней.
В этой работе предлагается именно такой путь аппроксимации моды распределения, а именно приближенное нахождение величин с, , на которых частичная сумма ряда Эджворта при учете центральных моментов ц4, ц5 и достигает максимального значения.
Далее приведем выражения для функций, аппроксимирующих плотность вероятности ф(£) при учете центрального моментов цр до порядка р включительно, подставляя выражения для полиномов Эрмита:
а) при р = 3:
ф£) «Ф(3)£) = ф(^)(Сз^3 - 3Сз^ +1);
б) при р = 4:
Ф(^) - ф(4) (£) = Ф(^)(Сб^6 + (С4 -15Сб + + С3 ^ +(45Сб-6С4);2-3С3 ^ + + 4С4 - 15Сб +1);
в) при р = 5:
ф(^) - ф(5)® = ф(^)(С9^9 + (С7 - 36С9)^7 +
+ Сб^6 + (С5 - 21С7 + 378С9)^5 + (С4 --15С6)^4 + (С3 - 10С5 + 105С7 - 1260С9)^3 + + (45С6 - 6С4)^2 + (15С5 - 3С3 - 105С7 +
+ 945С9^ + 3С4 -15Сб +1); г) при р = 6:
ф® - ф(бД) = Ф(^)(С!2^12 + (Сю - ббС^Г + + С£9 + (С8 - 45С10 + 1485С12)^8 + (С7 --36С9)^7 + (С6* - 28С8 + 630С10 -
- 13860С12)^6 + (С5 - 21С7 + 378С9 +
+ (С4 - 15С6* + 210С8 - 3150С10 + 51975С12 + + (С3 - 10С5 + 105С7 - 1260С9)^3 + (45С6* -
- 6С4 - 420С8 + 4725С10 - 62370С12)^2 +
+ (15С5 - 3С3 - 105С7 + 945С9 )£ + 3С4 - 15С6* + + 105С8 -945С10 + 10395С12 +1),
где
1 М3
С _-^3-, С _-1 -3
1 [ Мч
3! а
,3 ' 4
4! 1а4
Сб = -I , С = 1 -10М-
6 6! ^а3) 5 5! (а5 а3,
С _ 35 м,[м±-31 С _ 280ГМ3
4 7! а31а4 3), 4 9! (а3)
С*-1 -15 ^ + 30
С. I I —16 _4
6! V а
(2)
С8 _ 56 4 [М| -10 М|)+35 - 3 1 , 8 8! а3 и5 а3) 8! (а4 )
_ 2Ш0ГМзТ [ М4 - ^ С10 10! ^а31 (а4 31,
С _ 15400 Г^Т
С12 12! (а3) .
*
Для того чтобы найти величину , на которой частичная сумма ряда Эджворта достигает максимального значения, применим необходимые условия экстремума:
ё ф( р )(0
_ 0.
В результате получим алгебраическое уравнение, степень которого на единицу больше степени полинома, входящего в соответствующую частичную сумму:
а) при р = 3 это уравнение 4-й степени: Ф(3)(£) = -^4 + 6Сз^2 -^-3Сз = 0;
б) при р = 4 это уравнение 7-й степени:
Ф(4)® = -С6^7 + (21С6 -С4);5 -Сз^4 +
+ (10С4 - 105С6)^3 + 6С3+ + (105С6 - 15С4 -1)£- 3С3 = 0;
в) при р = 5 это уравнение 10-й степени:
Фо® = -С9Г + (45С9 - С7- С6^7 +
+ (28С7 -С5 - 630С9);6 + (21С6 - С4)^5 + + (15С5 -С3 -210С7 + 3150С9)^4 + + (10С4 - 105С6)^3 + (6С3 - 45С5 + + 420С7 - 4725С9 )^2 + (105С6 - 15С4 -1)£ + + 15С5 - 3С3 - 105С7 + 945С9 = 0;
г) при р = 6 это уравнение 13-й степени: Ф(6)® = -С^13 + (78С12 - С^11 - С9^10 +
+ (55С10 -С8 -2145С12)^9 + (45С9 -С7)^8 + + (36С8 -990С10 + 25740С12 - С*)^7 + + (28С7 -С5 - 630С9)^6 + (21С* - С4 --378С8 + 6930С10 - 135135С12)^5 + (15С5 -- С3 - 210С7 + 3150С9)^4 + (10С4 - 105С* + + 1260С8 - 17325С10 + 270270С12)^3 + + (6С3 - 45С5 + 420С7 - 4725С9)^2 + (105С* -- 15С4 -945С8 + 10395С10 - 135135С12 -1)£ + + 15С5 - 3С3 - 105С7 + 945С9 = 0.
Общее условие для всех приведенных выше алгебраических уравнений д3 Ф 0, т. е. ненулевой коэффициент асимметрии. Корни уравнения 4-й степени Ф(3)(£) = 0 можно отыскать, применяя известные соотношения [17], и в эквивалентной записи это уравнение применялось в [6, 7]. Для других уравнений предлагается применять метод Ла-герра [18, 19]. При решении модельных примеров применялась система компьютерной математики МаШсаё, а именно встроенная функция ро1утоо1& [20]. Для этой же цели могут быть использованы и другие методы, позволяющие эффективно находить корни алгебраических уравнений.
Даже для случая р = 6 количество корней алгебраического уравнения невелико и выбор того из них, на котором достигается
максимальное значение частичной суммы ряда Эджворта, можно осуществить с помощью простого перебора. Как только такой корень с, найден, аппроксимация моды распределения случайной величины X будет иметь вид
X * = М + а£*. (3)
При д3 = 0 нужно положить = 0 и X = М.
В качестве примеров рассмотрим некоторые типовые распределения случайных величин. Первые два из них - это показательное распределение и распределение х2 Они интересны тем, что формула (1) дает для них точное значение моды. Запишем плотность вероятности для показательного распределения с параметром X > 0:
ф(х) = , х > 0. Мода такого распределения равна нулю.
Плотность вероятности для распределения х с к степенями свободы определяется выражением
, , хкЛ-х/2 к
ф( х) =-, к = —, х >0,
^ 2кГ(к) 2
где Г(к) - гамма-функция. Мода равна к - 2 при к > 2.
Точные и приближенные значения моды для показательного распределения и распределения х , полученные по формулам (1) и (3), приведены в табл. 1 и 2. Точные значения выделены.
Таблица 1
Точные и приближенные значения моды для показательного распределения
X (1) (3) р = 3 (3) р = 4 (3) р = 5 (3) р = 6
0.5 0 0.985 0.655 0.493 0.394
1 0 0.493 0.327 0.247 0.197
2 0 0.246 0.164 0.123 0.099
5 0 0.099 0.065 0.049 0.039
10 0 0.049 0.033 0.025 0.020
Таблица 2
Точные и приближенные значения моды
2
для распределения х
к (1) (3) р = 3 (3) р = 4 (3) р = 5 (3) р = 6
2 0 0.985 0.655 0.493 0.394
3 1 1.850 1.485 1.313 1.212
4 2 2.753 2.373 2.206 2.117
5 3 3.678 3.295 3.140 3.067
10 8 8.463 8.116 8.029 8.007
20 18 18.289 18.033 18.004 18.000
Точность приближений моды распределения, полученных по формуле (3), невелика, особенно для показательного распределения, что можно объяснить тем, что соответствующая плотность вероятности имеет разрыв в нуле. Для распределения X формула (3) дает хорошее приближение при p = 6 и больших k, но это неудивительно, так как с ростом k это распределение приближается к нормальному.
Далее запишем плотность вероятности для распределения Рэлея с параметром о > 0 (этот параметр является модой распределения):
Ф( х) = "Г е
2а2
х > 0,
а также плотность вероятности для логарифмически нормального распределения с параметрами ц и о > 0 (мода распределения
1 (1п х -|1)2
Ф( х) =
равна ви
х > 0.
хс^Т/к
Приближенные значения моды для распределения Рэлея и логарифмически нормального распределения, полученные по формулам (1) и (3), приведены в табл. 3 и 4. Наиболее близкие к точному значения выделены.
Таблица 3 Приближенные значения моды для распределения Рэлея
о (1) (3) Р = 3 (3) р = 4 (3) Р = 5 (3) р = 6
0.5 0.523 0.538 0.507 0.504 0.502
1 1.047 1.076 1.013 1.007 1.004
2 2.093 2.153 2.026 2.014 2.008
4 4.186 4.306 4.053 4.029 4.017
8 8.373 8.611 8.105 8.057 8.033
Таблица 4 Приближенные значения моды для логарифмически нормального распределения
о (1) (3) Р = 3 (3) р = 4 (3) Р = 5 (3) р = 6
0.5 0.3 1.473 1.536 1.512 1.512 1.506
0.5 0.4 1.294 1.468 1.425 1.43 1.378
0.5 0.5 0.997 1.387 1.335 1.369 1.070
1 0.5 1.644 2.288 2.202 2.257 1.764
2 0.5 4.468 6.218 5.985 6.136 4.795
Здесь можно отметить, что точность приближений моды на основе формулы (3)
выше уже при p > 4 для распределения Рэлея. А для логарифмически нормального распределения наилучший результат достигается для формулы (3) при p = 4 или p = 6 (здесь формула (1) дает наименее точный результат оценивания).
Отметим, что с помощью частичных сумм ряда Эджворта можно построить и более точные приближения к плотности вероятности, следовательно, и к моде распределения. Однако это вряд ли целесообразно, поскольку противоречит основной идее: снижению вычислительных затрат. С увеличением максимального порядка учитываемых в разложении центральных моментов будет, очевидно, увеличиваться и время расчетов. В такой ситуации может оказаться, что приближенное нахождение моды по оценке плотности вероятности, например по гистограмме или ядерной оценке [14], более предпочтительно как по точности, так и по вычислительным затратам.
ЗАДАЧА ОПТИМАЛЬНОЙ ФИЛЬТРАЦИИ
Перейдем к задаче оптимальной фильтрации и рассмотрим модель объекта наблюдения, описываемую стохастическим дифференциальным уравнением Ито [3, 10, 11]:
dX ^) = f (^ X ^)) сЧ + , X ^)) dW ^),
X (О = X0,
в котором XеЯ" - вектор состояния; tеT = = t ] - отрезок времени функционирования системы; х): Т*Я"^Я", х):
- заданные функции; W(t) - 5-мерный стандартный винеровский процесс.
Модель измерительной системы записывается в форме
СУ (г) = с(г, X (г)) си + ^) СУ (г), (5)
где УеЯт - вектор измерений; с(^ х): Т*Я"
^Ят, С(0: Т^ЯтЫ - заданные функции
т
(для произведения С(0С (t) существует обратная матрица q(t)); ¥(?) - С-мерный стандартный винеровский процесс (№(?), У({) и X0 независимы).
Задача оптимальной фильтрации состоит
в нахождении оценки Хф по результатам измерений У^ ={У(т), те[^, t]} . При использовании МАП-критерия [3]
(4)
2
х
X(t) = arg max ф(/, х 1YJ"),
xeR"
где ф(t, х | Y0t) - апостериорная плотность вероятности вектора состояния X, т. е. оптимальная оценка - это мода апостериорного распределения.
Заметим, что мода апостериорного распределения, вообще говоря, может определяться неоднозначно, но здесь в качестве оптимальной оценки по МАП-критерию берется значение вектора состояния, на котором достигается глобальный максимум апостериорной плотности вероятности.
СТАТИСТИЧЕСКИЙ АЛГОРИТМ ФИЛЬТРАЦИИ
В одном из предыдущих разделов рассмотрены различные варианты приближенного нахождения моды распределения. Используем их, а именно формулу (3), как часть статистического алгоритма фильтрации, в основу которого положен алгоритм непрерывного фильтра частиц. Дополнительным упрощением служит нахождение координат оценки с помощью формулы (3) по моментам, которые соответствуют маргинальным плотностям вероятности [7].
В примерах, приведенных выше, рассматривались некоторые типовые распределения случайных величин. Однако в задаче фильтрации такие варианты для распределения если и возможны, то только для начального состояния. Апостериорное распределение даже для модельных задач фильтрации сложно идентифицировать, оно меняется с течением времени. Поэтому выбор конкретного варианта приближенного нахождения моды распределения нужно делать по результатам тестирования алгоритма фильтрации для конкретной задачи.
Далее приведем алгоритм приближенного решения задачи оптимальной фильтрации по МАП-критерию с применением ряда Эджворта. В основе алгоритма лежит непрерывный фильтр частиц, т. е. предполагается моделирование пары случайных процессов, первый из них X(t) задается уравнением (4), второй процесс ro(t) определяется следующим образом [4]:
ш(>) = ехр ст (т, X (т))д(тМ7 (т) --1Г 'с т(т, X (т))д(т)с(', X (т)М т}.
2 -"о )
Для реализации алгоритмов строятся дискретные аппроксимации случайных процессов Х(') и ю(') с помощью известных численных методов решения стохастических дифференциальных уравнений [4, 15].
Подобные статистические алгоритмы фильтрации приведены в работах [4, 13, 21], только в качестве критерия оптимальности использовался традиционный критерий минимума среднеквадратической ошибки оценивания [3]. В [7] алгоритм сформирован на основе МАП-критерия с нахождением оценки по формуле (1).
Шаг 1 (инициализация). Задать М - число моделируемых траекторий пары процессов (Х(У),ю(У)), т. е. число частиц; И - шаг численного интегрирования. Выбрать максимальный порядок р учитываемых центральных моментов: 3 <р < 6. Получить реализации векторов Х'0 согласно заданному
распределению начального вектора состояния Х0 (начальные положения частиц). Положить юг0 = 1, I = 1,2,...М и к = 0.
Шаг 2 (приближенное нахождение оптимальной оценки в текущем узле сетки по времени). Оценить апостериорные математические ожидания и апостериорные центральные моменты координат вектора состояния с учетом измерений , т. е. по реализациям {Хк 2 М с весовыми коэффициентами {шк}г=1,2,...,м :
1 M
M *, у = , у,
k i=1
П м I
6к,у !>k(Xk,у -Mk,у^
r, k, у
1M
= ^ (Хк,у - Хк, у )Г,
Ч i=1
где j = 1,2,...," - номер координаты, r = 3,...,p, fik - сумма всех весов:
м
Q k=s
®k.
i =1
Найти оценку моды апостериорного распределения, выбирая для ее каждой координаты корень уравнения Ф(р)(£) = 0, на котором соответствующая частичная сумма ряда Эджворта ф(р)(£) достигает максимального значения (числовые коэффициенты (2), задающие функции Ф(р>(^) и ф(р)(£), выражаются через оценки апостериорных центральных моментов). Результатом для каж-
е*
дого 7 является величина с)к 7, которая преобразуется по формуле (3) в 7-ю координату оценки моды апостериорного распределения в момент времени tk:
Ъ,1 = Мk, 7 * 7 Й,, .
Проверить условие Т - tk = 0. Если оно выполнено, то завершить процесс. Положить г = 1.
Шаг 3 (моделирование траекторий объекта наблюдения и соответствующих весовых коэффициентов). Получить реализацию вектора состояния в следующем узле сетки tk+l = tk + Ь и обновить весовой коэффициент:
X;+1 = ^ , X;, Ь), <+1 = < ехр {ст (tk, Xk)q(tk)(У- У(tk)) -
h
c \tk, ) q (tk )c (tk, Xk)}.
Шаг 4 (формирование циклов по времени и частицам). Если i = M, то положить k = = k + 1 и перейти к шагу 2; если i < M, то положить i = i + 1 и перейти к шагу 3.
В приведенном алгоритме предполагается, что шаг численного интегрирования выбран таким образом, что T/h - целое число. Функция F(t, X, h) определяется выбранным численным методом решения стохастических дифференциальных уравнений, например для стохастического метода Эйлера
F (t, X, h) = X + hf (t, X ) + 4h a(t, X )AW, где AW - 5-мерный случайный вектор, координаты которого независимы и имеют стандартное нормальное распределение (для моделирования приращений стандартного ви-неровского процесса). При каждом обращении к этой функции нужно моделировать новую реализацию AW.
Вектор Mk с координатами Mkj, который определяется на шаге 2 алгоритма, -
это приближение к оптимальной оценке вектора состояния по критерию минимума среднеквадратической ошибки оценивания.
На шаге 2 можно оценить начальные моменты:
r ,k, j
1 N
TT" К (Xk, j )r
"k i=1
при г = 2,..., р, выразив затем через них центральные моменты, используя известные соотношения [8].
ПРИМЕР РЕШЕНИЯ ЗАДАЧИ ОПТИМАЛЬНОЙ ФИЛЬТРАЦИИ
Пусть X(t) - скалярный случайный процесс, задаваемый простейшим уравнением: dX(г) = 0, X(0) = X,,,
где ^Т = [0,1], X0 - случайная величина, имеющая стандартное нормальное распределение.
Задача состоит в нахождении оценки Х^) по результатам измерений У при условии, что модель измерительной системы имеет вид СУ(г) = с(7 (t) - X(t))dt + СУ(г), У(0) = 0,
О и« и«
где с(х) = 24 + 6х + 3хт, 7($) = X + t, X -реализация случайной величины X0, У(t) -стандартный одномерный винеровский процесс.
Такая модель системы наблюдения и оценивания характерна для задач оценки погрешности навигационной системы по данным карты геофизического поля [1] (задача идентификации параметра X"). В частности, она использовалась в [7].
Результаты оценивания для значения параметра X = -1 приведены на рис. 1. При моделировании использовался стохастический метод Эйлера (для измерительной системы, для объекта наблюдения X(t) = X0), число частиц N = 102, Ь = 0.01 - шаг численного интегрирования. Дополнительно мода оценивалась по гистограмме, которая строилась на отрезке [-3,3] с шагом 0.06. Далее на рис. 2 и 3 показаны результаты оценивания для того же значения параметра X = -1, но при N = 103 и N = 104 соответственно.
Из приведенных результатов видно, что оценка моды по гистограмме чувствительна к выбору числа частиц М, когда оно невелико (для стабилизации результата требуется увеличение М).
№ - параметр X
* - оценка моды по гистограмме - оценка по асимметрии Шарлье ххх ~ п0 частичной сумме ряда Эджворта (р = 3) +_1_+ - по частичной сумме ряда Эджворта (р = 4) ВЭ€! - по частичной сумме ряда Эджворта (р = 5) 0 0 0 - по частичной сумме ряда Эджворта (р = 6)
Рис. 1. Результаты оценивания при М = 102
Оценка моды по центральным моментам апостериорного распределения в меньшей степени зависит от М, т. е. подтверждается утверждение о том, что оценка апостериорной плотности вероятности по сравнению с оценкой моментов требует увеличения числа моделируемых траекторий.
Рис. 2. Результаты оценивания при M = 10
Рис. 3. Результаты оценивания при M = 104
Стоит отметить, что с ростом максимального порядка p для учитываемых центральных моментов для частичной суммы ряда Эджворта в оценке моды все больше проявляется погрешность вычислений, что приводит к заметным осцилляциям. В рассматриваемом примере это заметно для случаев p = 5 и p = 6. В такой ситуации можно ограничиться значением p = 4.
На всех рисунках также приведены результаты приближенного нахождения моды апостериорного распределения по формуле (1), т. е. с помощью нахождения асимметрии Шарлье [7]. Дополнительно отметим, что оценка моды по гистограмме зависит не только от числа частиц M, но и от конкретных реализаций их начального положения, задаваемых на шаге 1 алгоритма. Центральные моменты апостериорного распределения также зависят от начального положения частиц, но в меньшей степени.
Все вычисления проведены с помощью специализированного пакета программ для системы компьютерной математики Mathcad [22].
ЗАКЛЮЧЕНИЕ
В работе на основе непрерывного фильтра частиц предлагается решение задачи оптимальной фильтрации по МАП-критерию.
С помощью фильтра частиц в предлагаемом алгоритме оцениваются апостериорные центральные моменты до 3-го, 4-го, 5-го или 6-го порядков включительно, а по ним приближенно оценивается мода апостериорного распределения. Для этого используется аппроксимация апостериорной плотности вероятности частичными суммами ряда Эджворта. Апробация предложенного алгоритма проведена на модельном примере идентификации параметра.
СПИСОК ЛИТЕРАТУРЫ
1. Степанов О. А. Основы теории оценивания с приложениями к задачам обработки навигационной информации: в 2 т. СПб.: ЦНИИ «Электроприбор», 2017. [ O. A. Stepanov, Fundamentals of estimation theory with applications to navigation information processing, (in Russian). Saint-Petersburg: TsNII "Elektropribor", 2017. ]
2. Bar-Shalom Y., Li X. R., Kirubarajan T. Estimation with applications to tracking and navigation. N. Y.: John Wiley & Sons, 2001. 581 p. [ Y. Bar-Shalom, X. R. Li, T. Kirubarajan, Estimation with applications to tracking and navigation. N. Y.: John Wiley & Sons, 2001. ]
3. Пантелеев А. В., Руденко Е. А., Бортаковский А. С. Нелинейные системы управления: описание, анализ и синтез. М.: Вузовская книга, 2008. 312 с. [ A. V. Panteleev, E. A. Rudenko, A. S. Bortakovskii, Nonlinear control systems: description, analysis and synthesis, (in Russian). Moscow: Vuzovskaya Kniga, 2008. ]
4. Bain A., Crisan D. Fundamentals of stochastic filtering. Springer, 2009. 394 p. [ A. Bain, D. Crisan, Fundamentals of stochastic filtering. Springer, 2009. ]
5. Косачев И. М., Чугай К. Н., Рыбаков К. А. Методология высокоточной нелинейной фильтрации случайных процессов в стохастических динамических системах с фиксированной структурой. Часть 1 // Труды МАИ. 2019. № 105. [ I. M. Kosachev, K. N. Chugai, K. A. Rybakov, "Methodology of high-precision nonlinear filtering of random processes in stochastic dynamical systems with a fixed structure (Part 1)," (in Russian), in Trudy MAI, no. 105, 2019. ]
6. Косачев И. М., Чугай К. Н., Рыбаков К. А. Методология высокоточной нелинейной фильтрации случайных процессов в стохастических динамических системах с фиксированной структурой. Часть 2 // Труды МАИ. 2019. № 106. [ I. M. Kosachev, K. N. Chugai, K. A. Rybakov, "Methodology of high-precision nonlinear filtering of random processes in stochastic dynamical systems with a fixed structure (Part 2)," (in Russian), in Trudy MAI, no. 106, 2019. ]
7. Chugai K., Kosachev I., Rybakov K. Approximate MMSE and MAP estimation using continuous-time particle filter // AIP Conference Proceedings. 2019. Vol. 2181. Id 020001. [ K. Chugai, I. Kosachev, K. Rybakov, "Approximate MMSE and MAP estimation using continuous-time particle filter", in AIP Conference Proceedings, vol. 2181, id 020001, 2019. ]
8. Cramer H. Mathematical methods of statistics. Princeton University Press, 1999. 575 p. [ H. Cramer, Mathematical methods of statistics. Princeton University Press, 1999. ]
9. Математическое обеспечение для анализа многомерных нелинейных стохастических систем / В. С. Пугачев
и др. // Автоматика и телемеханика. 1991. № 1. С. 87-97. [ V. S. Pugachev, et. al., "Software for the analysis of multidimensional nonlinear stochastic systems," (in Russian), in Avtomatika i telemekhanika, no. 1, pp. 87-97, 1991. ]
10. Синицын И. Н. Фильтры Калмана и Пугачева. М.: Логос, 2007. 776 с. [ I. N. Sinitsyn, Kalman and Pugachev filters, (in Russian). Moscow: Logos, 2007. ]
11. Косачев И. М., Ерошенков М. Г. Аналитическое моделирование стохастических систем. Минск: Наука и техника, 1993. 264 с. [ I. M. Kosachev, M. G. Eroshenkov, Analytical modeling of stochastic systems, (in Russian). Minsk: Nauka i tekhnika, 1993. ]
12. Косачев И. М. Методология высокоточной нелинейной фильтрации случайных процессов в стохастических динамических системах с фиксированной структурой // Вестник Воен. акад. Респ. Беларусь. 2014. № 4 (45). С. 125-161. [ I. M. Kosachev, "Methodology of high-precision nonlinear filtering of random processes in stochastic dynamical systems with a fixed structure," (in Russian), in Vestnik Voen. akad. Resp. Belarus', no. 4 (45), pp. 125-161, 2014. ]
13. Рыбаков К. А., Ющенко А. А. Непрерывные фильтры частиц и их реализация в реальном масштабе времени // Вестник ВГУ. Серия: Системный анализ и информационные технологии. 2018. № 3. С. 56-64. [ K. A. Rybakov, A. A. Yushchenko, "Continuous particle filters and its real-time implementation", (in Russian), in Vestnik VGU. Seriya: Sis-temnyi analiz i informatsionnye tekhnologii, no. 3, pp. 56-64, 2018. ]
14. Silverman B. W. Density estimation for statistics and data analysis. N. Y.: Chapman & Hall, 1986. 184 p. [ B. W. Silverman, Density estimation for statistics and data analysis. N. Y.: Chapman & Hall, 1986. ]
15. Аверина Т. А. Статистическое моделирование решений стохастических дифференциальных уравнений и систем со случайной структурой. Новосибирск: Изд-во СО РАН, 2019. 350 с. [ T. A. Averina, Statistical modeling of solutions of stochastic differential equations and systems with a random structure, (in Russian). Novosibirsk: Izd-vo SO RAN, 2019. ]
16. Млечин В. В. Теория радиоэлектронного преодоления. Анализ воздействия помех на радиотехнические системы и устройства. М.: Радиотехника, 2009. 976 с. [ V. V. Mlechin, Theory of radio-electronic overcoming. Analysis of the impact of interference on radio systems and devices, (in Russian). Moscow: Radiotekhnika, 2009. ]
17. Korn G. A., Korn T. M. Mathematical handbook for scientists and engineers. Dover Publ., 2000. 1152 p. [ G. A. Korn, T. M. Korn, Mathematical handbook for scientists and engineers. Dover Publ., 2000. ]
18. Pan V. Y. Solving a polynomial equation: some history and recent progress // SIAM Review. 1997. Vol. 39, no. 2. P. 187-220. [ V. Y. Pan, "Solving a polynomial equation: some history and recent progress," in SIAM Review, vol. 39, no. 2, pp. 187-220, 1997. ]
19. Мёллер Х. Алгоритм Лагерра/сумм степеней для эффективной и надежной аппроксимации всех корней многочлена // Проблемы передачи информации. 2015. Т. 51, № 4. С. 60-70. [ H. Moller, "The Laguerre-and-sums-of-powers algorithm for the efficient and reliable approximation of all polynomial roots,"(in Russian), in Problemy peredachi informacii, vol. 51, no. 4, pp. 60-70, 2015. ]
20. Очков В. Ф. Mathcad 14 для студентов и инженеров: русская версия. СПб.: БХВ-Петербург, 2009. 512 с.
[ V. F. Ochkov, Mathcad 14 for students and engineers: Russian version, (in Russian). Saint-Petersburg: BHV-Petersburg, 2009. ]
21. Рыбаков К. А. Статистические алгоритмы оптимальной фильтрации сигналов в нелинейных диффузионно-скачкообразных стохастических системах // Вестник УГАТУ. 2016. Т. 20, № 4 (74). С. 107-113. [ K. A. Rybakov, "Statistical algorithms of optimal filtering problem for nonlinear jump-diffusion models," (in Russian), in Vestnik UGATU, vol. 20, no. 4 (74), pp. 107-113, 2016. ]
22. Кудрявцева И. А., Руденко Е. А., Рыбаков К. А. Программное обеспечение оптимального оценивания состояний стохастических динамических систем // Информационные и телекоммуникационные технологии. 2019. № 43. С. 23-28. [ I. A. Kudryavtseva, E. A. Rudenko, K. A. Rybakov, "Software for optimal state estimation in stochastic dynamical systems," (in Russian), in Informatsionnyye i telekommu-nikatsionnyye tekhnologii, no. 43, pp. 23-28, 2019. ]
ОБ АВТОРЕ
РЫБАКОВ Константин Александрович, доц. каф. математической кибернетики. Дипл. математик-инж. (МАИ, 2002). Канд. физ.-мат. наук по сист. анализу, управл. и обраб. информ. (МАИ, 2006). Иссл. в обл. методов моделирования, анализа и синтеза стохастических систем управления.
METADATA
Title: Approximate solution of optimal filtering problem with respect to the MAP estimation using partial sums of the Edgeworth series. Author: K. A. Rybakov Affiliation:
Moscow Aviation Institute (National Research University) (MAI), Russia. Email: rkoffice@mail.ru Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 24, no. 1 (87), pp. 93-102, 2020. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: Based on the continuous-time particle filter and approximation of the probability density by partial sums of the Edgeworth series, a statistical algorithm for the approximate solution of the optimal filtering problem with respect to the MAP estimation is formed. In the proposed algorithm, the mode of the conditional distribution for each unobservable state vector component is expressed in terms of conditional third, fourth, fifth or sixth order central moments. Numerical example is given. Key words: maximum of conditional probability density; mode of distribution; Hermite polynomials; optimal estimation; optimal filtering; Edgeworth series; statistical algorithm; stochastic differential equation; stochastic system; particle filter.
About author:
RYBAKOV, Konstantin Alexandrovich, Associate Prof., Dept. of Math. Cybernetics. Dipl. Math.-Engineer (MAI, 2002). Cand. of Phys.-Math. Sci. (MAI, 2006).
