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

Многочасти́чный фильтр^[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].

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

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

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

x_{k}=f_{k}(x_{k-1},v_{k})

,

где $f_{k}$ функция изменения состояния системы, $v_{k}$ — случайная величина, возмущающее воздействие.

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

y_{k}=h_{k}(x_{k},w_{k})

,

где $h_{k}$ функция измерения, $w_{k}$ — случайная величина, шум измерений.

Функции $f_{k}$ и $h_{k}$ в общем случае нелинейные, а статистические характеристики шума системы ( $v_{k}$ ) и измерений ( $w_{k}$ ) предполагаются известными.

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

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

Рассмотрим дискретный марковский процесс $\{X_{n}\}_{n\geqslant 1}$ со следующими распределениями вероятностей:

X_{1}\sim \mu (x_{1})\quad

и

X_{n}\mid (X_{n-1}=x_{n-1})\sim f(x_{n}\mid x_{n-1})

,

(1)

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

Здесь нотация $X\mid Y\sim f(\dots )$ означает, что $X$ при условии $Y$ распределено как $f(\dots )$ .

Реализации процесса $\{X_{n}\}$ (скрытые переменные $x_{n}$ ) наблюдаются посредством другого случайного процесса $\{Y_{n}\}_{n\geqslant 1}$ — процесса измерений — с маргинальными плотностями:

Y_{n}\mid (X_{n}=x_{n})\sim h(y_{n}\mid x_{n})

,

(2)

где $h(y_{n}\mid x_{n})$ — условная плотность вероятности (плотность измерений), измерения считаются статистически независимыми.

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

{\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}}

Для простоты считаем, что переходная плотность и плотность измерений не зависят от $n$ . Параметры модели считаются заданными.

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

Уравнение (1) определяет априорное распределение для процесса $\{X_{n}\}$ :

p(x_{1:n})=\mu (x_{1})\prod _{k=2}^{n}f(x_{k}\mid x_{k-1})

(3)

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

p(y_{1:n}\mid x_{1:n})=\prod _{k=1}^{n}h(y_{k}\mid x_{k})

,

(4)

Здесь и далее нотация $x_{k:l}$ для $k\leqslant l$ обозначает $(x_{k},\dots ,x_{l})$ .

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

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})}}

,

(5)

где (здесь $dx_{1:n}$ — доминирующая мера):

p(y_{1:n})=\int p(x_{1:n})p(y_{1:n}\mid x_{1:n})\,dx_{1:n}

.

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

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

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

{\bar {\theta }}=\int \theta (x)p(x)\,dx

,

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

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

{\bar {\theta }}\approx {\frac {1}{N}}\sum _{i=1}^{N}\theta (x^{(i)})

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

w_{i}={\frac {r(x^{(i)})}{\sum _{j=1}^{N}r(x^{(j)})}}

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

{\bar {\theta }}=\int \theta (x)r(x)q(x)\,dx\approx \sum _{i=1}^{N}w_{i}\theta (x^{(i)})

,

Перевыборка

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

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

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

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

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

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

f_{n}(x_{1:n})={\frac {\phi _{n}(x_{1:n})}{Z_{n}}}

, где

\phi _{n}\colon {\mathcal {X}}^{n}\to \mathbb {R} ^{+}

известна поточечно, а

Z_{n}=\int \phi _{n}(x_{1:n})\,dx_{1:n}

— нормализующая, возможно неизвестная, постоянная, то

SMC-алгоритм будет находить приближения $f_{k}(x_{1:k})$ и оценки $Z_{k}$ для $k=1,2,\dots$ .

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

\phi _{n}(x_{1:n})=p(x_{1:n})p(y_{1:n}\mid x_{1:n})

и

Z_{n}=p(y_{1:n})

,

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

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})

.

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

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})

— предиктор,

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})}}

— корректор.

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

Алгоритм

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

   Алгоритм МЧФ
   -- инициализация
   для i = 1...N:
     выборка  $\xi _{0}^{(i)}$  из  $q_{0}(x_{0}\mid y_{0})$ 
     -- начальные веса
      $\omega _{0}^{(i)}:=h(y_{0}\mid \xi _{0}^{(i)})\mu (\xi _{0}^{(i)})\ /\ q_{0}(\xi _{0}^{(i)}\mid y_{0})$  
   кц
   для n = 1...T:
     если ПЕРЕВЫБОРКА то
       -- выбор индексов  $j_{i}\in \{1,\dots ,N\}$  N частиц в соответствии с весами
        $j_{1:N}$  = SelectByWeight( $\{w_{n-1}^{(j)}\}$ )
       для i = 1...N:
          $x_{n-1}^{(i)}:=\xi _{n-1}^{(j_{i})}$ 
          $w_{n-1}^{(i)}:=1/N$ 
     иначе
       для i = 1...N:
          $x_{n-1}^{(i)}:=\xi _{n-1}^{(i)}$ 
     для i = 1...N:
       -- шаг распространения частицы
        $\xi _{n}^{(i)}\sim q_{n}(\xi _{n}^{(i)}\mid \xi _{n-1}^{(i)},y_{n})$ 
       -- обновление весов
        $\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})$  
     кц
     -- нормализация весов
      $s:=\sum _{j=1}^{N}\omega _{n}^{(j)}$ 
     для i = 1...N:
        $w_{n}^{(i)}:=\omega _{n}^{(i)}/s$ 
   кц

См. также

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

Примечания

↑ ¹ ² Микаэльян, 2011.
↑ Gordon, Salmond, Smith, 1993.
↑ ¹ ² ³ ⁴ ⁵ ⁶ ⁷ ⁸ 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. Архивировано 4 марта 2016 года.
↑ ¹ ² ³ 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] ¹ ² Микаэльян, 2011.

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

[_b007ae1902e5143a-3] ¹ ² ³ ⁴ ⁵ ⁶ ⁷ ⁸ 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. Архивировано 4 марта 2016 года.

[_3265dc2a6a5a4019-5] ¹ ² ³ 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]