Многочастичный фильтр

Многочасти́чный фильтр^[1] (МЧФ, англ. particle filter — «фильтр частиц», «частичный фильтр», «корпускулярный фильтр») — последовательный метод Монте-Карло — рекурсивный алгоритм для численного решения проблем оценивания (фильтрации, сглаживания), особенно для нелинейных и не-гауссовских случаев. Со времени описания в 1993 году^[2] Н. Гордоном, Д. Салмондом и А. Смитом используется в различных областях — навигации, робототехнике, компьютерном зрении.

В сравнении с обычно применяемыми для подобных задач методами — расширенными фильтрами Кальмана (EKF) — многочастичные фильтры не зависят от методов линеаризации или апроксимации. Обычный EKF плохо справляется с существенно нелинейными моделями, а также в случае шумов системы и измерений, сильно отличающихся от гауссовых, поэтому были разработаны различные модификации, такие как UKF (англ. unscented KF), QKF (англ. Quadrature KF) и т. п.^[3]. Следует отметить, что в свою очередь многочастичные фильтры более требовательны к вычислительным ресурсам.

Термин «particle filter» был дан Дел Моралом в 1996 году^[4], а «sequential Monte Carlo» — Лю (Liu) и Ченом (Chen) в 1998.

Многие используемые на практике многочастичные фильтры выводятся применением последовательного метода Монте-Карло к последовательности целевых распределений^[5].

Постановка задачи

МЧФ предназначен для оценки последовательности скрытых переменных [math]\displaystyle{ x_n }[/math] для [math]\displaystyle{ n = 1, 2, \dots }[/math] на основании наблюдений [math]\displaystyle{ y_n }[/math] при [math]\displaystyle{ n = 1, 2, \dots }[/math]. Для простоты изложения будем считать, что рассматривается динамическая система, и [math]\displaystyle{ x_n }[/math] и [math]\displaystyle{ y_n }[/math] — действительные вектора состояния и измерений соответственно^[1].

Стохастическое уравнение состояния системы имеет вид:

[math]\displaystyle{ x_k = f_k(x_{k-1}, v_k) }[/math],

где [math]\displaystyle{ f_k }[/math] функция изменения состояния системы, [math]\displaystyle{ v_k }[/math] — случайная величина, возмущающее воздействие.

Уравнение измерений:

[math]\displaystyle{ y_k = h_k(x_k, w_k) }[/math],

где [math]\displaystyle{ h_k }[/math] функция измерения, [math]\displaystyle{ w_k }[/math] — случайная величина, шум измерений.

Функции [math]\displaystyle{ f_k }[/math] и [math]\displaystyle{ h_k }[/math] в общем случае нелинейные, а статистические характеристики шума системы ([math]\displaystyle{ v_k }[/math]) и измерений ([math]\displaystyle{ w_k }[/math]) предполагаются известными.

Задачей фильтрации является получение оценки [math]\displaystyle{ \hat{x}_k }[/math] на основе известных к моменту [math]\displaystyle{ k }[/math] результатов измерений [math]\displaystyle{ y_{1:k} }[/math].

Скрытая марковская модель и байесовский вывод

Рассмотрим дискретный марковский процесс [math]\displaystyle{ \{X_n\}_{n \geqslant 1} }[/math] со следующими распределениями вероятностей:

[math]\displaystyle{ X_1 \sim \mu(x_1)\quad }[/math] и
[math]\displaystyle{ X_n \mid (X_{n-1} = x_{n-1}) \sim f(x_n \mid x_{n-1}) }[/math],

(1)

где [math]\displaystyle{ \mu(x) }[/math] — плотность вероятности, [math]\displaystyle{ f(x_n \mid x_{n-1}) }[/math] — условная плотность вероятности (переходная плотность вероятности) при переходе от [math]\displaystyle{ x_{n-1} }[/math] к [math]\displaystyle{ x_n }[/math].

Здесь нотация [math]\displaystyle{ X \mid Y \sim f(\dots) }[/math] означает, что [math]\displaystyle{ X }[/math] при условии [math]\displaystyle{ Y }[/math] распределено как [math]\displaystyle{ f(\dots) }[/math].

Реализации процесса [math]\displaystyle{ \{X_n\} }[/math] (скрытые переменные [math]\displaystyle{ x_n }[/math]) наблюдаются посредством другого случайного процесса [math]\displaystyle{ \{Y_n\}_{n \geqslant 1} }[/math] — процесса измерений — с маргинальными плотностями:

[math]\displaystyle{ Y_n \mid (X_n = x_n) \sim h(y_n \mid x_n) }[/math],

(2)

где [math]\displaystyle{ h(y_n \mid x_n) }[/math] — условная плотность вероятности (плотность измерений), измерения считаются статистически независимыми.

Модель может проиллюстрирована следующей диаграммой переходов:

[math]\displaystyle{ \begin{array}{cccccccccc} X_1&\rightarrow&X_2&\rightarrow&X_3&\rightarrow&X_4&\rightarrow&\ldots&\\ \downarrow&&\downarrow&&\downarrow&&\downarrow&&\ldots&\\ Y_1&&Y_2&&Y_3&&Y_4&&\ldots& \end{array} }[/math]

Для простоты считаем, что переходная плотность и плотность измерений не зависят от [math]\displaystyle{ n }[/math]. Параметры модели считаются заданными.

Определённая таким образом модель системы и измерений известна как скрытая марковская модель^[6].

Уравнение (1) определяет априорное распределение для процесса [math]\displaystyle{ \{X_n\} }[/math]:

[math]\displaystyle{ p(x_{1:n}) = \mu(x_1) \prod^n_{k=2} f(x_k \mid x_{k-1}) }[/math]

(3)

Аналогично (2) задаёт функцию правдоподобия:

[math]\displaystyle{ p(y_{1:n}\mid x_{1:n}) = \prod^n_{k=1} h(y_k\mid x_k) }[/math],

(4)

Здесь и далее нотация [math]\displaystyle{ x_{k:l} }[/math] для [math]\displaystyle{ k \leqslant l }[/math] обозначает [math]\displaystyle{ (x_k, \dots, x_l) }[/math].

Таким образом, байесовский вывод для [math]\displaystyle{ \{X_{1:n}\} }[/math] при известных реализациях измерений [math]\displaystyle{ \{Y_{1:n}\} }[/math], обозначенных соответственно как [math]\displaystyle{ \{x_{1:n}\} }[/math] и [math]\displaystyle{ \{y_{1:n}\} }[/math], будет опираться на апостериорное распределение

[math]\displaystyle{ p(x_{1:n}\mid y_{1:n}) = \frac{p(x_{1:n}) p(y_{1:n}\mid x_{1:n})}{p(y_{1:n})} }[/math],

(5)

где (здесь [math]\displaystyle{ dx_{1:n} }[/math] — доминирующая мера):

[math]\displaystyle{ p(y_{1:n}) = \int p(x_{1:n}) p(y_{1:n}\mid x_{1:n})\,dx_{1:n} }[/math].

Выборка по значимости

См. также Выборка по значимости.

Метод Монте-Карло позволяет оценивать свойства довольно сложных распределений вероятностей, например, путём вычисления средних и дисперсии в виде интеграла^[3]:

[math]\displaystyle{ \bar{\theta} = \int \theta(x) p(x)\,dx }[/math],

где [math]\displaystyle{ \theta(x) }[/math] — функция для оценивания. Например, для среднего можно положить: [math]\displaystyle{ \theta(x) = x }[/math].

В случае невозможности аналитического решения, задача может быть решена численно генерированием случайных выборок с плотностью [math]\displaystyle{ p(x) }[/math], обозначим их как [math]\displaystyle{ {x^{(i)}}_{1\leqslant i \leqslant N} }[/math], и получением среднего арифметического по точкам выборки^[3]:

[math]\displaystyle{ \bar{\theta} \approx \frac{1}{N} \sum_{i=1}^N \theta(x^{(i)}) }[/math]

В более общем случае, когда выборка из [math]\displaystyle{ p }[/math] затруднена, применяется другое распределение [math]\displaystyle{ q }[/math] (так называемое англ. instrumental or importance distribution), а для сохранения несмещённости оценки вводятся весовые коэффициенты [math]\displaystyle{ w_i }[/math] на основе отношения [math]\displaystyle{ r(x^{(i)}) = p(x^{(i)}) / q(x^{(i)}) }[/math]^[3]:

[math]\displaystyle{ w_i = \frac{r(x^{(i)})} {\sum^{N}_{j=1} r(x^{(j)})} }[/math]

после чего вычисляет взвешенное среднее:

[math]\displaystyle{ \bar{\theta} = \int \theta(x) r(x) q(x)\,dx \approx \sum_{i=1}^N w_i \theta(x^{(i)}) }[/math],

Перевыборка

Хотя вспомогательное распределение используется в основном для упрощения выборки из основного распределения [math]\displaystyle{ p }[/math], часто применяется процедура «выборки и перевыборки по значимости» (англ. sampling importance resampling, SIR). Эта процедура состоит из двух этапов: собственно выборки по значимости с вычислением весов [math]\displaystyle{ w_i }[/math], и дополнительной выборки точек, учитывающих эти веса^[3].

Перевыборка особенно необходима для последовательных фильтров^[3].

Последовательный метод Монте-Карло

Методы многочастичной фильтрации и сглаживания являются наиболее известными примерами алгоритмов последовательного метода Монте-Карло (англ. sequential Monte Carlo, SMC). До такой степени, что в литературе часто не делают между ними различия. Тем не менее, SMC включает в себя более широкий класс алгоритмов, применимых для описания более сложных приблизительных методов фильтрации и сглаживания^[7].

Последовательные методы Монте-Карло являются классом методов Монте-Карло, которые производят последовательную выборку из последовательности целевых плотностей вероятностей [math]\displaystyle{ \{f_n(x_{1:n})\} }[/math] увеличивающейся размерности, где каждое [math]\displaystyle{ f_n(x_{1:n}) }[/math] определено на декартовой степени [math]\displaystyle{ \mathcal{X}^n }[/math]^[5].

Если записать плотность как:^[5]

[math]\displaystyle{ f_n(x_{1:n}) = \frac{\phi_n(x_{1:n})}{Z_n} }[/math], где

[math]\displaystyle{ \phi_n\colon \mathcal{X}^n \to \mathbb{R}^{+} }[/math] известна поточечно, а

[math]\displaystyle{ Z_n = \int \phi_n(x_{1:n})\,dx_{1:n} }[/math] — нормализующая, возможно неизвестная, постоянная, то

SMC-алгоритм будет находить приближения [math]\displaystyle{ f_k(x_{1:k}) }[/math] и оценки [math]\displaystyle{ Z_k }[/math] для [math]\displaystyle{ k = 1, 2, \dots }[/math].

Например, для случая фильтрации можно положить (см. (5)):

[math]\displaystyle{ \phi_n(x_{1:n}) = p(x_{1:n}) p(y_{1:n}\mid x_{1:n}) }[/math] и

[math]\displaystyle{ Z_n = p(y_{1:n}) }[/math],

из чего будем иметь:

[math]\displaystyle{ f_n(x_{1:n}) = \frac{p(x_{1:n}) p(y_{1:n}\mid x_{1:n})}{p(y_{1:n})} = p(x_{1:n}|y_{1:n}) }[/math].

Опуская вывод, схему предиктор-корректор можно представить в следующем виде^[3]:

[math]\displaystyle{ p(x_{1:n}\mid y_{1:n-1}) = p(x_{1:n-1}\mid y_{1:n-1}) f(x_n\mid x_{n-1}) }[/math] — предиктор,

[math]\displaystyle{ p(x_{1:n}\mid y_{1:n}) = \frac{h(y_n\mid x_n) p(x_{1:n}\mid y_{1:n-1})}{p(y_n\mid y_{1:n-1})} }[/math] — корректор.

Множитель [math]\displaystyle{ (p(y_n\mid y_{1:n-1}))^{-1} }[/math] — нормализующая постоянная, которая не требуется для обычного SMC-алгоритма.

Алгоритм

Типичный алгоритм многочастичного фильтра можно представить в следующем виде^[3]:

   Алгоритм МЧФ
   -- инициализация
   для i = 1...N:
     выборка [math]\displaystyle{ \xi_0^{(i)} }[/math] из [math]\displaystyle{ q_0(x_0\mid y_0) }[/math]
     -- начальные веса
     [math]\displaystyle{ \omega_0^{(i)} := h(y_0\mid \xi_0^{(i)}) \mu(\xi_0^{(i)})\ /\ q_0(\xi_0^{(i)}\mid y_0) }[/math] 
   кц
   для n = 1...T:
     если ПЕРЕВЫБОРКА то
       -- выбор индексов [math]\displaystyle{ j_i \in \{1,\dots,N\} }[/math] N частиц в соответствии с весами
       [math]\displaystyle{ j_{1:N} }[/math] = SelectByWeight([math]\displaystyle{ \{ w_{n-1}^{(j)}\} }[/math])
       для i = 1...N:
         [math]\displaystyle{ x_{n-1}^{(i)} := \xi_{n-1}^{(j_i)} }[/math]
         [math]\displaystyle{ w_{n-1}^{(i)} := 1 / N }[/math]
     иначе
       для i = 1...N:
         [math]\displaystyle{ x_{n-1}^{(i)} := \xi_{n-1}^{(i)} }[/math]
     для i = 1...N:
       -- шаг распространения частицы
       [math]\displaystyle{ \xi_{n}^{(i)} \sim q_n(\xi_{n}^{(i)}\mid \xi_{n-1}^{(i)}, y_n) }[/math]
       -- обновление весов
       [math]\displaystyle{ \omega_n^{(i)} := w_{n-1}^{(i)} h(y_n\mid \xi_n^{(i)}) f(\xi_n^{(i)}\mid x_{n-1}^{(i)})\ /\ q_n(\xi_n^{(i)}\mid x_{n-1}^{(i)}, y_n) }[/math] 
     кц
     -- нормализация весов
     [math]\displaystyle{ s := \sum^{N}_{j=1} \omega_{n}^{(j)} }[/math]
     для i = 1...N:
       [math]\displaystyle{ w_{n}^{(i)} := \omega_{n}^{(i)} / s }[/math]
   кц

См. также

Фильтр Кальмана#UKF

Примечания

↑ ^1,0 ^1,1 Микаэльян, 2011.
↑ Gordon, Salmond, Smith, 1993.
↑ ^3,0 ^3,1 ^3,2 ^3,3 ^3,4 ^3,5 ^3,6 ^3,7 Cappé, Godsill, Moulines, 2007.
↑ Del Moral, Pierre. Non Linear Filtering: Interacting Particle Solution. (англ.) // Markov Processes and Related Fields. — 1996. — Vol. 2, no. 4. — P. 555–580.
↑ ^5,0 ^5,1 ^5,2 Doucet, Johansen, 2011.
↑ Doucet, Johansen, 2011, 2.1 Hidden Markov Models and Inference Aims.
↑ Doucet, Johansen, 2011, 3 Sequential Monte Carlo Methods.

Литература

Doucet Arnaud, Johansen Adam M. A Tutorial on Particle Filtering and Smoothing: Fifteen Years Later // The Oxford Handbook of Nonlinear Filtering / D. Crisan, B. Rozovsky. — Oxford : Oxford University Press, 2011. — P. 656—704. — ISBN 978-0-19-953290-2.
Cappé, Olivier and Godsill, Simon J. and Moulines, Eric. An Overview of Existing Methods and Recent Advances in Sequential Monte Carlo // Proceedings of the IEEE. — IEEE, 2007. — Т. 95, № 5. — P. 899—924. — ISSN 0018-9219. Архивировано 10 марта 2016 года.

Doucet, Arnaud and de Freitas, Nando and Gordon, Neil. An Introduction to Sequential Monte Carlo Methods // Sequential Monte Carlo Methods in Practice / Doucet, Arnaud and de Freitas, Nando and Gordon, Neil. — Springer New York. — 3-14 p. — ISBN 978-1-4419-2887-0.
Arulampalam, M.S. and Maskell, S. and Gordon, N. and Clapp, T. A Tutorial on Particle Filters for Online Nonlinear/non-Gaussian Bayesian Tracking (англ.) // Trans. Sig. Proc.. — IEEE Press, 2002. — Vol. 50, no. 2. — P. 174—188. — ISSN 1053-587X. См. также более раннюю версию (англ.)
Gordon, N.J.; Salmond, D.J.; Smith, A.F.M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation (англ.) // IEEE Proceedings F, Radar and Signal Processing. — IET, 1993. — Vol. 140, no. 2. — P. 107—113. — doi:10.1049/ip-f-2.1993.0015.
Микаэльян С. В. Методы фильтрации на основе многоточечной аппроксимации плотности вероятности оценки в задаче определения параметров движения цели при помощи измерителя с нелинейной характеристикой // Наука и образование : электронное издание. — МГТУ им. Н. Э. Баумана, 2011. — ISSN 1994-0408. Архивировано 4 марта 2016 года.
Ristic, B., Arulampalam, S., Gordon, N. Beyond the Kalman Filter — Particle Filters for Tracking Applications. — Artech House, 2004. — 299 p. — ISBN 9781580536318.

Simon, Dan. 15 The particle filter // Optimal State Estimation: Kalman, H_∞, and Nonlinear Approaches. — Wiley-Interscience, 2006. — P. 461—480. — ISBN 0471708585.

Ссылки

Particle Filter, SciPy Cookbook

[_1332e6b1235634ae-1] 1,0 ^1,1 Микаэльян, 2011.

[_09e8c8b3705a5607-2] Gordon, Salmond, Smith, 1993.

[_b007ae1902e5143a-3] 3,0 ^3,1 ^3,2 ^3,3 ^3,4 ^3,5 ^3,6 ^3,7 Cappé, Godsill, Moulines, 2007.

[4] Del Moral, Pierre. Non Linear Filtering: Interacting Particle Solution. (англ.) // Markov Processes and Related Fields. — 1996. — Vol. 2, no. 4. — P. 555–580.

[_3265dc2a6a5a4019-5] 5,0 ^5,1 ^5,2 Doucet, Johansen, 2011.

[_5821243940019ece-6] Doucet, Johansen, 2011, 2.1 Hidden Markov Models and Inference Aims.

[_6ae5e8bdcaccbef9-7] Doucet, Johansen, 2011, 3 Sequential Monte Carlo Methods.

[1]

[2]

[3]

[4]

[5]

[6]

[7]