ИЗВЕСТИЯ РАН. ТЕОРИЯ И СИСТЕМЫ УПРАВЛЕНИЯ, 2013, № 3, с. 16-27
УПРАВЛЕНИЕ В СТОХАСТИЧЕСКИХ СИСТЕМАХ И В УСЛОВИЯХ НЕОПРЕДЕЛЕННОСТИ
УДК 519.245 519.226.3
ИССЛЕДОВАНИЕ ПОГРЕШНОСТИ ВЫЧИСЛЕНИЯ ОПТИМАЛЬНОЙ БАЙЕСОВСКОЙ ОЦЕНКИ МЕТОДОМ МОНТЕ-КАРЛО В НЕЛИНЕЙНЫХ ЗАДАЧАХ* © 2013 г. Н. А. Берковский, О. А. Степанов
Санкт-Петербург, ЦНИИ "Электроприбор", Санкт-Петербургский государственный политехнический ун-т Поступила в редакцию 04.06.12 г., после доработки 28.09.12 г.
Предложена процедура оценивания погрешности вычисления оптимальной байесовской оценки, основанная на исследовании вида функции плотности распределения вероятности для отношения двух интегралов, вычисляемых с помощью метода Монте-Карло. Выявлена связь предложенного метода с известным дельта-методом. Проведены результаты численных экспериментов для тестовой задачи и задачи навигации по точечным ориентирам, подтвердившие эффективность предложенного метода и правильность полученных выводов.
DOI: 10.7868/S0002338813010034
Введение. При обработке навигационной информации широко используются методы оптимального байесовского оценивания [1—7]. Хорошо известно, что оптимальная байесовская оценка с минимальной дисперсией представляет собой математическое ожидание, соответствующее апостериорной плотности для оцениваемых параметров [5, 6]. Простые алгоритмы вычисления оптимальной оценки в виде рекуррентных соотношений калмановского типа могут быть получены только в линейных гауссовских задачах. Однако на практике, в том числе и при обработке навигационной информации, нередко приходится сталкиваться с нелинейными задачами. В этом случае для нахождения оптимальных оценок для решения навигационных задач в реальном времени на борту подвижных объектов необходимы специальные упрощенные (субоптимальные) алгоритмы. На сегодняшний день разработан целый арсенал различных методов построения таких алгоритмов. Исследования в этой области продолжаются, в том числе и в связи с увеличением потребности решения нелинейных задач обработки навигационной информации применительно к различным подвижным объектам нового класса [7—14]. Наиболее корректная оценка эффективности алгоритмов в рамках байесовского подхода основана на сопоставлении точности, достигаемой с их помощью, с точностью оптимального алгоритма. Эта точность определяется, как правило, с помощью определения выборочных значений среднеквадратической ошибки оптимальной оценки путем многократного расчета оптимальной оценки. При нахождении этой оценки в камеральных условиях объем вычислений играет второстепенную роль. Существенным здесь является умение определять факт достижения заданной точности. Известно, что оптимальная байесовская оценка при решении нелинейных задач может быть представлена в виде частного двух интегралов [5, 6]. Для их вычисления широкое применение получил метод Монте-Карло и его различные модификации [5, 6, 15—17]. Одно из достоинств метода Монте-Карло заключается в возможности анализировать точность в ходе вычисления. При этом хорошо известно, как контролировать погрешность приближенного вычисления одного интеграла-числителя или знаменателя [18]. Для оценки же точности вычисления частного на практике обычно используют приближенный метод, основанный на гауссовской аппроксимации частного и получивший наименование в литературе дельта-метода [19, 20].
В настоящей работе рассматривается более общий подход, базирующийся на аппроксимации функции плотности распределения частного в предположении, что числитель и знаменатель имеют приближенно гауссовское распределение с соответствующими параметрами. Обоснована законность использования такой аппроксимации. Разработана вычислительная схема, применяющая предложенную аппроксимацию для контроля погрешности расчета байесовской оценки.
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-08-00372a).
Приведены результаты численных экспериментов. Сопоставлены результаты, полученные при помощи предложенного алгоритма, дельта-метода и простейшей формулы для оценки погрешности вычисления частного.
1. Постановка задачи. Пусть известна совместная функция плотности распределения вероятности (далее — плотность) /(х, у) вектора состояния х е Я" и вектора измерений у е Ят. В этом случае оптимальная байесовская оценка для x может быть представлена в виде [5, 6]
| xf (x, y )dx
x(y) =
R"
| f (x, y )dx
Я"
Для 1-й компоненты х' вектора состояния эту оценку можно записать как
| xf (x, y )dx
_ _ R"_
1(y) Jf (x,y)dx '
(y) = ^ = ^-, i = 1,..., П. (1.1)
Я"
Будем полагать, что №е приближение к байесовской оценке (1.1) вычисляется с помощью метода Монте-Карло
-4 (у) = Л (у) /Л (у),
где
~1 { ) = 1 у4£М 12(у) = 1 уАхь!), Л ) Р(хк) , 2 ™ ^ Р(хк) ,
а хь х2, ...хх — последовательность независимых случайных векторов в Я", распределенных с некоторой плотностью р^), неотрицательной на носителе функции / (х, у), рассмотренной при фиксированном у. Это соответствует известной схеме существенной выборки [18] для получения статистических оценок интегралов числителя и знаменателя в (1.1). В численных примерах вид р(^) конкретизируется, но основные результаты статьи, приведенные в разд. 2 и 3, не зависят от
выборар^). При этом хк = (х\,х2к,..., х"). В дальнейшем зависимость от у в записи оценок 1Х (у)
и 12 (у) иногда для краткости опускается. Согласно дельта-методу в качестве аппроксимации для
распределения статистики Х'н используется нормальное распределение с плотностью
N
x, X' (y), X' (y)2
v_2 _ 2 2cov(ibi2
+ ^2 _
V Л 12 I1I2 JJ
(1.2)
где ст1 и ст2 — среднеквадратические отклонения 11 и I2 [19, 20]. Здесь и далее N(x, a, b) — плотность нормального распределения с аргументом x (скалярным или векторным), математическим ожиданием a и ковариационной матрицей b, а X ~ N (a, b) означает, что случайная величина или вектор X имеет нормальное распределение с математическим ожиданием a и ковариационной матрицей b, при этом в скалярном случае ковариационная матрица b — это просто дисперсия X. Векторный аргумент выделяется жирным шрифтом.
2. Аппроксимация распределения оценки частного двух интегралов при ее вычислении методом Монте-Карло. Сначала рассмотрим следующую вспомогательную задачу. Пусть имеются две случайные величины 11 и I2 такие, что распределение двумерного случайного вектора I = (I2)
гауссовское, т.е. I ~ (1,2), где I = (1,12)T и
Z =
2
01 РСТ1СТ 2
VPCT1CT2 J
Здесь 11 и 12 — математические ожидания 11 и 12; , а2, р — среднеквадратические отклонения и коэффициент корреляции для 11 и 12. Рассмотрим случайную величину & = 11 / 12.
Утверждение 1. Пусть случайный вектор (12) ~ N(I, Е). Тогда случайная величина & = 11/12 будет распределена с плотностью g (^), определяемой по формуле
8 ^ =ЖпеХР ^ 2 ^ " 121)2■ 2 2,1 ^ (*)' (2
где
2(а2 - 2ра!а2г + а2г2),
а
(г)
-(I2at +1г - pgigi (( + IiZ)) [j - 2F ( ^а? +1^г - ра^ (i + IiZ) ^
(а? - 2pCTiCT2г + а2г2)1,5
V"1U2
2 т . /т 2 т чч2Л
aiaW(1 - p2)(а2 - 2pgiCT2г + а22г2)
+ + (2.2)
2аа ./i expI (12^ - Лра1а2 + г(11^2 - I2pCTiCT2)) 2а1а 2V1 -P expl —- 2 2---
¿(г) =_у 2(1 -рУ^а2 - 2pgiCT2г + 2) у. (2.3)
(а2 - 2ра1а2г + а2г2>/2П
Здесь F (г) = -у== | exp j dt — функция Лапласа.
Доказательство утверждения 1 приведено в Приложении.
Отметим без доказательства некоторые свойства плотности g(z).
1. В случае I2 Ф 0, ст2 > 0 |р| < 1 (что выполнено во всех ситуациях, рассматриваемых в статье) функция g(z) непрерывна и на бесконечности g (г) = 0(г ~2). Непрерывность следует из положительности выражения (а2 - 2ра1а2г + а2г2) = (ст1 - ра2г)2 + (1 - р2)^2г2, входящего во все знаменатели соотношений (2.1)-(2.3). Оценка на бесконечности следует из того, что неотрицательные выражения под экспонентами в (2.1) и (2.3) при г ^ да стремятся к конечным пределам, соответственно, экспоненты не стремятся к нулю и убывание g(z) на бесконечности определяется первым множителем в (2.2) и знаменателем в (2.3). Таким образом, вероятностное распределение с плотностью g(z) является распределением с "тяжелыми хвостами". Следовательно, у случайных величин, распределенных с плотностью g(z), математическое ожидание и дисперсия в смысле интеграла Стильтьеса не существуют [21].
2. Если а1 = да2, где q > 0 — некоторый коэффициент, то при г Ф I1 /I2 и фиксированных значениях параметров I1, I2 Ф 0 и |р| < 1 справедливо lim g(г) = 0. Это происходит в силу стремле-
ст2 ^ 0
ния к нулю экспонент в (2.1) и (2.3), которые убывают быстрее, чем возрастает первый множитель в (2.2). Соотношение а1 = да 2 будем считать выполненным во всех рассматриваемых в статье приложениях функции g(z), так как q — не что иное, как квадратный корень из отношения дисперсий статистических оценок числителя и знаменателя по методу Монте-Карло.
3. При а1 = да 2 и достаточно малых значениях а2 функция g(z) имеет колоколообразную форму с ярко выраженным максимумом в точке г max, расположенной в окрестности точки г = I1 /I2. При ст2 ^ 0 имеет место сходимость гmax ^ I1/12. На рис. 1 показаны графики функции g(z) при различных значениях а2, при ст1 = ст2 (д = 1) и при I1 = 2, I2 = 1, р = 0.5. Из рис. 1 хорошо видно, что при уменьшении величины а2 (а вместе с ней и а1) график функции g(z) становится симметричным, расположение максимума приближается к I1 /I2 = 2, значение экстремума увеличивается, а область существенно отличных значений от нуля уменьшается.
4. Допустим, что I1/I2 > 0. Тогда если значение коэффициента q близко к величине I1/I2, то при уменьшении а 2 имеет место увеличение максимума функции g(z) по мере того, как величина р приближается к единице. Если же значение q значительно отличается от I1 /I2, то форма графика g(z) слабо зависит от величины р. Это иллюстрируют графики на рис. 2, где функция g(z) рассмотрена как функция двух переменных z и р. Графики получены при фиксированных значениях I1 = 2, I2 = 1, ст2 = 0.1, при этом коэффициент д = 7 (рис. 2, а) и д = 2.5 (рис. 2, б). Хорошо
г
4
3 2 1 0 -1
! ! I 1 1 V а2 = 0. 1 05 !
\а2 = 0.1
а2 = 0.5
1 1 1 1 1 1 1
0.5 1.0 1.5 2.0 2.5 3.0 3.5
г
Рис. 1. Графики функции g (г) при трех различны
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.