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

Многочасти́чный фильтр^[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 ⩾ 1}$ со следующими распределениями вероятностей:

X_{1} \sim μ (x_{1})

и

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

,

(1)

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

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

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

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

,

(2)

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

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

\begin{array}{cccccccccc} X_{1} & \to & X_{2} & \to & X_{3} & \to & X_{4} & \to & \dots \\ ↓ & ↓ & ↓ & ↓ & \dots \\ Y_{1} & Y_{2} & Y_{3} & Y_{4} & \dots \end{array}

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

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

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

p (x_{1 : n}) = μ (x_{1}) \prod_{k = 2}^{n} f (x_{k} ∣ x_{k - 1})

(3)

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

p (y_{1 : n} ∣ x_{1 : n}) = \prod_{k = 1}^{n} h (y_{k} ∣ x_{k})

,

(4)

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

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

p (x_{1 : n} ∣ y_{1 : n}) = \frac{p (x_{1 : n}) p (y_{1 : n} ∣ x_{1 : n})}{p (y_{1 : n})}

,

(5)

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

p (y_{1 : n}) = \int p (x_{1 : n}) p (y_{1 : n} ∣ x_{1 : n}) d x_{1 : n}

.

Выборка по значимости[править]

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

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

\bar{θ} = \int θ (x) p (x) d x

,

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

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

\bar{θ} \approx \frac{1}{N} \sum_{i = 1}^{N} θ (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{θ} = \int θ (x) r (x) q (x) d x \approx \sum_{i = 1}^{N} w_{i} θ (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})$ определено на декартовой степени $𝒳^{n}$ ^[5].

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

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

, где

ϕ_{n} : 𝒳^{n} \to ℝ^{+}

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

Z_{n} = \int ϕ_{n} (x_{1 : n}) d x_{1 : n}

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

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

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

ϕ_{n} (x_{1 : n}) = p (x_{1 : n}) p (y_{1 : n} ∣ x_{1 : n})

и

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

,

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

f_{n} (x_{1 : n}) = \frac{p (x_{1 : n}) p (y_{1 : n} ∣ x_{1 : n})}{p (y_{1 : n})} = p (x_{1 : n} | y_{1 : n})

.

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

p (x_{1 : n} ∣ y_{1 : n - 1}) = p (x_{1 : n - 1} ∣ y_{1 : n - 1}) f (x_{n} ∣ x_{n - 1})

— предиктор,

p (x_{1 : n} ∣ y_{1 : n}) = \frac{h (y_{n} ∣ x_{n}) p (x_{1 : n} ∣ y_{1 : n - 1})}{p (y_{n} ∣ y_{1 : n - 1})}

— корректор.

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

Алгоритм[править]

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

   Алгоритм МЧФ
   -- инициализация
   для i = 1...N:
     выборка  $ξ_{0}^{(i)}$  из  $q_{0} (x_{0} ∣ y_{0})$ 
     -- начальные веса
      $ω_{0}^{(i)} : = h (y_{0} ∣ ξ_{0}^{(i)}) μ (ξ_{0}^{(i)}) / q_{0} (ξ_{0}^{(i)} ∣ 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)} : = ξ_{n - 1}^{(j_{i})}$ 
          $w_{n - 1}^{(i)} : = 1 / N$ 
     иначе
       для i = 1...N:
          $x_{n - 1}^{(i)} : = ξ_{n - 1}^{(i)}$ 
     для i = 1...N:
       -- шаг распространения частицы
        $ξ_{n}^{(i)} \sim q_{n} (ξ_{n}^{(i)} ∣ ξ_{n - 1}^{(i)}, y_{n})$ 
       -- обновление весов
        $ω_{n}^{(i)} : = w_{n - 1}^{(i)} h (y_{n} ∣ ξ_{n}^{(i)}) f (ξ_{n}^{(i)} ∣ x_{n - 1}^{(i)}) / q_{n} (ξ_{n}^{(i)} ∣ x_{n - 1}^{(i)}, y_{n})$  
     кц
     -- нормализация весов
      $s : = \sum_{j = 1}^{N} ω_{n}^{(j)}$ 
     для i = 1...N:
        $w_{n}^{(i)} : = ω_{n}^{(i)} / s$ 
   кц

См. также[править]

Фильтр Кальмана#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]