Kalman-and-Bayesian-Filters-in-Python 2．離散ベイズフィルター(1)

2.1 離散ベイズフィルター

今回はベイズフィルターについてです。
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/02-Discrete-Bayes.ipynb
実はカルマンフィルターはベイズフィルターの一種。
ベイズの概念のような確率論がどのようにフィルタリング&トラッキングに利用されているのかをお話ししていきます。

犬を追う～prior:事前確率～

notebookでは著者の好みにより犬のトラッキングを題材にするようです。ワンちゃんのいる可能性がある場所が仮に10か所で、それぞれ1/10の確率で存在するなら

In（事前確率）

import numpy as np
belief = np.array([1./10]*10)
print(belief)
belief

Out（事前確率）

[ 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1]

こんなかんじでnumpyを使って表現します。

廊下で～categorical distribution～

なぜか急に題材が変わり、三つのドア（センサー付き）がある廊下でおじさんがどこにいるのかという問題を定式化すると

from kf_book.book_plots import figsize, set_figsize
import kf_book.book_plots as book_plots
import matplotlib.pyplot as plt

#ドアの並び
belief = np.array([1./3, 1./3, 0, 0, 0, 0, 0, 0, 1/3, 0])
plt.figure()
set_figsize(y=2)
book_plots.bar_plot(belief)

おじさんの存在する確率分布

ドアの前にいると、必ずセンサーが反応すると仮定して、ある時刻にセンサーが反応するならば上図のような確率分布で表現できます。

この分布の様子を、categorical distribution(n回の観測による離散的な確率)と呼びます。

センサーから情報を抽出する

もしおじさんのセンサーが
・door
・move right
・door
という値を返したら、おじさんの居場所を推測することは可能でしょうか？
ーはい、可能です。
そのような値を返せるのは2番目の場所のみですから。
このようにして観測値を通して実際の状況を推定するというのが最終的に私たちがしたいことです。

でも現実は甘くない

極端な話ですが、1/3センサーが正しく動いてくれるとは限りません。なので他の場所に1/10の確率でいることを加味すると

beliefの更新

def update_belief(hall, belief, z, correct_scale):
    for i, val in enumerate(hall):
        if val == z:
            belief[i] *= correct_scale

belief = np.array([0.1] * 10)
reading = 1 # 1 is 'door'
update_belief(hallway, belief, z=reading, correct_scale=3.)
print('belief:', belief)
print('sum =', sum(belief))
plt.figure()
book_plots.bar_plot(belief)

Out

belief: [ 0.3  0.3  0.1  0.1  0.1  0.1  0.1  0.1  0.3  0.1]
sum = 1.6

確率の合計を１にするため正規化すると

normalization

belief / sum(belief)

Out

array([0.188, 0.188, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.188,
       0.062])

ただ、”3回正しくもあり誤りでもある”が併存している状態があるのはおかしいので”センサーが正しく動作する確率”も考慮してあげます

scale =  \frac{prob_{correct}}{prob_{incorrect}} = \frac{prob_{correct}} {1-prob_{correct}}

一般に、pythonで正規化する場合は以下の関数を使うことができます

nomalization

from filterpy.discrete_bayes import normalize
normalize(belief)

で、numpyを使うときにforループはなるべく使いたくないので、ブール代数を使ってドアの判定をしてあげます

hallway == 1

array([ True,  True, False, False, False, False, False, False,  True,
       False])

あらためてbeliefの更新をかけると

beliefの更新Ⅱ

from filterpy.discrete_bayes import normalize

def scaled_update(hall, belief, z, z_prob): 
    scale = z_prob / (1. - z_prob)
    belief[hall==z] *= scale
    normalize(belief)

belief = np.array([0.1] * 10)
scaled_update(hallway, belief, z=1, z_prob=.75)

print('sum =', sum(belief))
print('probability of door =', belief[0])
print('probability of wall =', belief[2])
book_plots.bar_plot(belief, ylim=(0, .3))

sum = 1.0
probability of door = 0.1875
probability of wall = 0.06249999999999999

あらかじめ与えられた
・belief（事前確率)
・hall（観測）
の情報を統合してあらたな状態を表現することが出来ました
こうして得られた結果をposterior:事後確率といいます

ちなみに、

belief[hall==z] *= scale

こちらの式で計算されているのはlikelihood:尤度と呼ばれるものです。(確率ではないので、1以上になりうる。何もなければ1をかけてそのまんまに）尤度を用いて事後確率を表現するなら

\mathtt{posterior} = \frac{\mathtt{likelihood} \times \mathtt{prior}}{\mathtt{normalization}}

あらためて更新式を書き直すと

bleiefの更新式Ⅲ

def scaled_update(hall, belief, z, z_prob): 
    scale = z_prob / (1. - z_prob)
    likelihood = np.ones(len(hall)) #np.one(すべての要素が1の行列を生成）
    likelihood[hall==z] *= scale
    return normalize(likelihood * belief)

・likelihood:尤度
・prior:事前確率
から
・posterior:事後確率
をreturnしていますね

参考

こちらのjupyter notebookをダウンロードして使用しています。

Kalman-and-Bayesian-Filters-in-Python
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python