はじめに
逐次ベイズフィルタの基本的な概要は極めて単純です。しかし非常に強力です。
制御の分野では遥か昔から状態観測器としてカルマンフィルタとして知られる逐次ベイズフィルタが有効活用されてきました。また、数理モデルによる演繹的なシミュレーションと、観測データによる機能的な推測を統合したデータ同化と呼ばれる分野でも、主にパーティクルフィルタが強力なツールとして利用されています。また自己位置推定、SLAMなど近年の自律移動ロボット技術に欠かせない物となっています。
必要なパーツ
逐次ベイズフィルタに必要なパーツは下記の通り、たったの２つです。これらを紹介する前に記法について整理しておきましょう。
時刻 $t$ での状態を $x _ t$ とし、観測を $y _ t$ とします。また、$x _ {1: t}$ と書いた場合には $1, ... , t$ までの $x$ すべて、すなわち $x _ 1, x _ 2, ..., x _ t$ の略記とします。$p(a)$ は $a$ の確率密度関数で、$p(a \mid b)$ は $b$ を所与としたときの $a$ の条件付き確率密度関数です。さて、では、逐次ベイズフィルタに必要な２つのパーツを見てみましょう。
予測
$$ p(x _ t \mid y _ {1 : t - 1}) $$
これは時刻 $t-1$ までの観測値を所与としたときの（観測データから見て）次の状態 $x _ t$ の確率密度です。これを使うことで、１つ先の未来の状態を確率的に表現できます。
観測更新
$$ p(x _ t \mid y _ {1 : t}) $$
これは時刻 $t$ までの観測値を所与としたときの（観測データから見て）現在の状態 $x _ t$ の確率密度です。これを使うことで、得られた観測から今まさにどういう状態であるのかを確率的に表現できます。
逐次ベイズフィルタの流れ
前提
さて、２つのパーツを使って具体的に行うのは、予測と観測更新を交互に実行していくことだけです。一度、頭を整理しましょう。今は予測と観測更新が具体的にどのようなものかは分からなくても良いです。兎にも角にも、$p(x _ t \mid y _ {1 : t - 1})$ という密度関数と $p(x _ t \mid y _ {1 : t })$ という密度関数を、今"手元に持ち合わせている"という状況です。
そういったときに、観測が逐次的に得られていく状況を想定します。
流れ
まずは興味を持っている状態 $x$ は、初期状態 $x _ 1$ から始まります。これは既知かもしれないですし、未知であるかもしれません。ここでは一般的な態度を取り未知であることにします。しかし、初期状態がどのような値を持つのかを確率的には見積もれる、すなわち $p(x _ 1)$ は目処付ができることにします。
すると初期状態は
$$ p(x _ 1) $$
で推測されることになります。ここに、状態 $x _ 1$ に関連のある $y _ 1$ という値を観測することになります。そこで状態 $x _ 1$ がどのようになっているかを推測するには
$$ p(x _ 1 \mid y _ 1) $$
という条件付き確率密度を利用します。これは $t=1$ のときの観測更新の密度関数 $p(x _ t \mid y _ {1 : t})$ であり、これは所与なのですから、$x _ 1$ がどんな状態であるのかを確率的に見積もることができます。単に $p(x _ 1)$ と目処付していたところ、観測値 $y _ 1$ によって推測を更新しました。
さて、今 $y _ 1$ という観測値を得たことによって、次の時刻の状態 $x _ 2$ の推測を $p(x _ 2 \mid y _ 1)$ によって実施できるようになりました。なぜなら予測のための密度関数 $p(x _ t \mid y _ {1 : t - 1})$ が所与であり、今は $t = 2$ というケースに当てはまるからです。今、$x _ 2$ に関して
$$ p(x _ 2 \mid y _ 1) $$
という推測を立てることができるようになったのです。 そして、再び次の観測値 $y _ 2$ が得られるとどうでしょう。観測更新の密度関数を $t = 2$ のケースで利用して
$$ p(x _ 2 \mid y _ {1: 2}) $$
によって $x _ 2$ を推測することができます。また、予測の密度関数によって一歩先の $x _ 3$ を推測することもできるようになります。
$$ p(x _ 3 \mid y _ {1 : 2}) $$
まとめると、時刻 $t$ において私達のできる推測は、$y _ {t}$ を観測することによって状態 $x _ t$ を$p(x _ t \mid y _ {1 : t})$ によって推測することと、更にこれまでデータを考慮して次の状態 $x _ {t+1}$ を $p(x _ {t+1} \mid y _ {1 : t})$ によって推測することです。あるいは見方を少し変えると、時刻 $t$ での状態 $x _ t$ に関して、$y _ t$ が観測できていない状況下において推測を予測の密度関数
$$ p(x _ t \mid y _ {1 : t-1}) $$
で推測することになり、観測が得られた時点で
$$ p(x _ t \mid y _ {1 : t}) $$
によって推測を行うというふうに見ることもできます。"予測"と"観測更新"という言葉にフィットするのはこちらの見方かもしれませんね。そして刻を進めて $t+1$ での $x _ {t+1}$ に興味が移るのならば、観測 $y _ {t + 1}$ がまだ手に入っていないときは $p(x _ {t+1} \mid y _ {1 : t})$ を使い、手に入った段階で $p(x _ {t+1} \mid y _ {1 : t+1})$ を使う…と続いていくのです。
予測の密度関数をどう使うのか
さて、ここまでの話で少し疑問に思われた方もいるかもしれません。
「あれ、予測の方いらなくない？」と。
$x _ t$ に興味がある場合にその推測をしなければならないのなら、観測 $y _ t$ を得てから行えばよいのではないでしょうか。ここまでの説明では、予測の方の密度関数は時刻 $t - 1$ と $t$ の間、観測 $y _ t$ がまだ手に入っていないときのその場しのぎ的に使っているようにしか見えないからです。
この"予測の方はいらないのではないか"という疑問は、２つの例により否定されます。
まずは１つめ、実際の時刻に照らし合わせると $t-1$ と $t$ のラグがとても長いにもかかわらず、状態 $x$ に応じた意思決定を常に迫れるというケースがあります。今 $t - 1$ と $t$ の間は実際の時刻で言うところの1秒間であるとします。では状態 $x _ t$ はそんな1秒毎にしか変化しないのでしょうか。そうではありません。時刻の離散化は通常、人間側の都合によって行われます。
例えば、ここでセンサが搭載された自動運転車を想像しましょう。1秒間というのは自動運転車のカメラの動作周期であるとしましょう（普通、そんな遅くないけど）。さて、今、カメラに写っている世界の状態（通行人や信号、自分自身の車速など）は１秒に１回しか更新されないのですが、実際の世界はカメラが更新されるその間にも変化し続けています。仮に自動運転車の操作が0.1秒毎に司令可能ならば、$t - 1$ と $t$ の間に相当する$t - 0.9, t- 0.8, ..., t-0.1$ という時刻における状態をどのように見積もり、どのような操作指令をすべきでしょうか。
まさか、$t-1$ 時点での観測更新の密度関数 $p( x _ {t - 1} \mid y _ {t - 1})$ に従って状態を見積もり、その後の状態を一定であると考えるなんてありえないはずです。カメラの更新が来る間にも、刻一刻と状態を予測しながら操作を実施するはずです（１秒毎にしか目を開けられない運転手の気持ちになれば、目をつぶっている間、自身の車速と目を開けたときに写った通行人の歩行状態から、状態を推測し続けるはずです）。
すなわち、このケースでは次の観測が来るまでの時刻では常に予測の確率密度関数ひたすら使い続けることになります（予測と観測更新を交互にやると言いながら、若干発展的な例になってしまいましたが、制御の分野ではかなりオーソドックスな例です）。予測の密度関数は、応用上多くの場合必須なのです。
そして、もう１つの理由が、下記に述べるように、"観測更新の式に予測の式が必要になる"ことと"予測の式に観測更新の式が必要になる"という、ある種本質的な理由になります（じゃあ最初からそれを見せてくれ）。すなわち、そもそも両方が支え合っているのです。より具体的には状態 $x _ t$ に対する予測の式は
$$ p(x _ t \mid y _ {1 : t-1}) = f (p(x _ {t-1} \mid y _ {1 : t-1})) $$
といった具合に、$x _ {t-1}$ に対する観測更新の式を引用する形で記述されます。また、$x _ t$ に対する観測更新の式は
$$ p(x _ t \mid y _ {1 : t}) = g(p(x _ t \mid y _ {1 : t-1})) $$
と、$x _ t$ に対する予測の式を引用する形になります。なので、すでに説明したとおり、時刻が進むごとに両方を交互に実行する必要があるということです。その詳細は次に見ていきましょう。
各パーツの式展開
予測の密度関数
まず、周辺化公式 $p(a) = \int p(a, b) db$ というものを利用します。この周辺化公式は、別の確率変数 $c$ がいて $p(a \mid c)$ というものに対しても利用が可能で
$$ p(a \mid c) = \int p(a, b \mid c) db $$
と書くことができます。これは右から左に読むと分かりやすく、"確率変数 $b$ については積分して消しちゃってるから左辺では消えているよ"ということです。
これを用いると、予測の密度関数は
$$ p(x _ t \mid y _ {1 : t - 1}) = \int p(x _ t, x _ {t-1} \mid y _ {1: t - 1}) d x _ {t-1} $$
と書くことができます。さて、この積分の中身を更に見ていきましょう。
$$ p(x _ t, x _ {t-1} \mid y _ {1: t - 1}) = \frac{p(x _ t, x _ {t-1}, y _ {1: t - 1})}{p(y _ {1: t - 1}) } $$
と同時分布と条件付き分布の関係式 $p(a, b, c) = p(a, b \mid c)p(c)$ を使って変形ができます。この関係式の肝は、"条件側に移動した変数がいるなら、移動した変数だけの分布がくっついてくる”ということです。条件側に移動する変数は複数個でもよく$p(a, b, c) = p(a \mid b, c)p(b, c)$ などとしても良いのです。このことを利用して、上記右辺の分子を
$$ p(x _ t, x _ {t-1}, y _ {1: t - 1})=p(x _ t \mid x _ {t-1}, y _ {1: t-1})p(x _ {t-1}, y _ {1:t-1}) $$
と変形してしまいましょう。時刻 $t$ とそれ以前で分けるという具合です。更に右辺の第二因子に関しても同時分布と条件付き分布の関係式を使って状態と観測を分離してやることにします。
$$ p(x _ {t-1}, y _ {1:t-1}) = p(x _ {t-1} \mid y _ {1:t-1})p(y _ {1:t-1}) $$
と分解するととても都合が良いです。結局の所、同時分布に関しては
$$ p(x _ t, x _ {t-1}, y _ {1: t - 1})=p(x _ t \mid x _ {t-1}, y _ {1: t-1})p(x _ {t-1} \mid y _ {1:t-1})p(y _ {1:t-1}) $$
と分解されました。いま予測の密度関数の積分の中身が
$$ \begin{align} p(x _ t, x _ {t-1} \mid y _ {1: t - 1}) &= \frac{p(x _ t, x _ {t-1}, y _ {1: t - 1})}{p(y _ {1: t - 1}) } \\ &= \frac{p(x _ t \mid x _ {t-1}, y _ {1: t-1})p(x _ {t-1} \mid y _ {1:t-1})p(y _ {1:t-1})}{p(y _ {1: t - 1}) } \\ &= p(x _ t \mid x _ {t-1}, y _ {1: t-1})p(x _ {t-1} \mid y _ {1:t-1}) \\ & = p(x _ t \mid x _ {t-1})p(x _ {t-1} \mid y _ {1:t-1}) \end{align} $$
という形になりました。ここで、最後の式変形は、状態の時間発展 $x _ {t-1} \rightarrow x _ t $ は、我々がこれまでに観測した $y _ {1:t-1}$ に依存するわけがない、という仮定に基づくものです。単に時間発展している状態を外から観測しているだけであって、観測しようがしまいが、現在の状態は前の状態に基づいて決まる、と考えているということです（もしも状態が、観測の仕方、計測の仕方によって影響を受け、次の状態が計測そのものに依存するような事があるならば、最後の式変形はできないというわけだ。しかし今回はそんなケースは想定しないことにする）。
これを予測の密度関数の積分の中身に形式的に放り込んでやれば予測の密度関数の式が得られるわけですが、ここで重要な事実を確認しておきましょう。上記の式には $p(x _ {t-1} \mid y _ {1:t-1})$ という部分がくっついています。これは、時刻 $t$ から見れば、１個前の時刻の観測更新の密度関数になっています。これは頭の中で、
$$ x _ {t - 1} \sim p(x _ {t-1} \mid y _ {1:t-1}) $$
と１個前の状態を想像しつつ、その状態を使って
$$ x _ t \sim p(x _ t \mid x _ {t-1}) $$
現在の状態を推測するようなイメージになります。しかしこのような推測を単発で行うのは少々頼りない気がします。実際の予測の密度関数は積分の形をしており、
$$ p(x _ t \mid y _ {1:t-1}) = \int p(x _ t \mid x _ {t-1})p(x _ {t-1} \mid y _ {1:t-1}) d x _ {t - 1} $$
となるので、$x _ {t - 1} \sim p(x _ {t-1} \mid y _ {1:t-1})$ という１個前の状態の生成に関して、考えうる$x _ {t-1}$ に関する値を無限に細かく考慮しているということになります。”これはパラメータを $x _ {t-1}$ とし、所与のデータを $y _ {1:t} $ とし、そして予測対象を $x _ t$ としたときのベイズ予測分布そのものです"。
更新の密度関数
さて、状態 $x _ t$ に対する観測更新の密度関数を見ていきましょう。まずは下記のように観測データを現在の値とそれまでの値に区別して書いてみます。
$$ \begin{align} p(x _ t \mid y _ {1 : t}) &= \frac{p(x _ t, y _ {1: t})}{p(y _ {1:t})}\\ &= \frac{p(x _ t, y _ t, y _ {1: t-1})}{p(y _ {t}, y _ {1: t-1})} \\ &= \frac{p(x _ t, y _ t, y _ {1: t-1})}{p(y _ {t} \mid y _ {1: t-1})p(y _ {1: t-1})} \end{align} $$
まずは例のごとく同時分布と条件付き分布の関係性を使って、分子に同時分布を持ってきました。その後、現在の観測値 $y _ t$ を強調して略記していた $y _ {1: t}$ から取り出し、更に分母に関しては現在時刻と過去の値で分離しています。
次は分子の同時分布について更に詳しく見ていきましょう。まず、時の流れを模擬すると同時分布に出てくる変数のうち、$y _ t$ という観測値が得られるのが最後になるはずです。この確率変数が決まっていく流れを記述してやることにすると、同時分布は
$$ \begin{align} p(x _ t, y _ t, y _ {1: t-1}) &= p(y _ t \mid x _ {t} , y _ {1: t - 1})p(x _ {t} , y _ {1: t-1}) \\ &= p(y _ t \mid x _ {t} , y _ {1: t - 1})p(x _ {t} \mid y _ {1: t-1})p(y _ {1: t-1}) \\ &= p(y _ t \mid x _ {t})p(x _ {t} \mid y _ {1: t-1})p(y _ {1: t-1}) \end{align} $$
と変形できます。ここで最後の式変形は $x _ t$ と観測 $y _ {1:t-1}$ が両方所与の時、$y _ t$ の確率を見積もるに当たって必要なのは $x _ t$ だけであるという仮定を使いました。これは、$x _ t$ を何らかの方法で観測して $y _ t$ の値を得ているという、今想定している内容に合致しています。一方で、仮に、計測方法に何らかの問題があり、$y _ t$ の値が現在の状態 $x _ t$ のみならず、これまでの $y _ {1 : t - 1}$ の観測結果を引きずってしまう（例えば１つ前の過去で計測精度が非常に高い値を得ると、今回は１回休み的に値がバグる、みたいなことがあるならば、過去の観測値は今回の観測値を見積もるための情報になる。今回はそういうケースは扱っていない。）
さて、出来上がった同時分布を観測更新の式に乗せてやれば
$$
\begin{align}
p(x _ t \mid y _ {1 : t})
&= \frac{p(x _ t, y _ t, y _ {1: t-1})}{p(y _ {t} \mid y _ {1: t-1})p(y _ {1: t-1})} \\
&= \frac{p(y _ t \mid x _ {t})p(x _ {t} \mid y _ {1: t-1})p(y _ {1: t-1})}{p(y _ {t} \mid y _ {1: t-1})p(y _ {1: t-1})} \\
&= \frac{p(y _ t \mid x _ {t})p(x _ {t} \mid y _ {1: t-1})}{p(y _ {t} \mid y _ {1: t-1})}
\end{align}
$$
となります。分母は $x _ t$ には関係ない変数しか出てきておらず、$x _ t$ から見たら単に所与の値 $y _ t, y _ {1: t} $ だけで構成されており定数に過ぎません。したがって、さしあたって分布の形状を考える上では
$$ p(x _ t \mid y _ {1 : t}) \propto p(y _ t \mid x _ {t})p(x _ {t} \mid y _ {1: t-1}) $$
と表記しておきましょう。さて、重要なのは $p(x _ {t} \mid y _ {1: t-1})$ が $x _ t$ に対する予測の密度関数になっているということです。観測更新というのは、過去の観測 $y _ {1: t-1}$ だけでひとまず $x _ t$ を予測するということをしておいて、それに加え、 $x _ t $ で件付けられた $y _ t$ の密度関数、すなわち現在の観測の密度関数を乗ずることで得られます。"観測更新の密度関数は、予測の密度関数を事前分布に、観測を尤度関数にした事後分布と見なせます"。
まとめ
まとめると、
- 逐次ベイズフィルタは予測の密度関数と観測更新の密度関数の２つのパーツから成る。
- 時刻 $t$ での予測の密度関数は $t-1$ での観測更新の密度関数を使って計算できる。
- 時刻 $t$ での観測更新の密度関数は $t$ での予測の密度関数を使って計算できる。
- 予測の密度関数は観測更新の密度関数を事後分布としたベイズ予測分布である。
- 観測更新の密度関数は予測の密度関数を事前分布とした事後分布である。。
すなわち、適当な初期分布を使ってベイズ予測分布を構築し、観測が得られたときにはベイズ予測分布を事前分布に据えた事後分布の算出を行い、次はその事後分布によってベイズ予測分布を構築し、…
と繰り返すのが逐次ベイズフィルタというわけです。 予測の式と観測更新の式をもう一度確認しましょう。
予測は
$$ p(x _ t \mid y _ {1:t-1}) = \int p(x _ t \mid x _ {t-1})p(x _ {t-1} \mid y _ {1:t-1}) d x _ {t - 1} $$
となり、観測更新は
$$ p(x _ t \mid y _ {1 : t}) \propto p(y _ t \mid x _ {t})p(x _ {t} \mid y _ {1: t-1}) $$
となります。これらの式が所与なら逐次ベイズフィルタを実行することで観測 $ y $ から状態 $x$ を推測できるのでした。
すると、応用上必要になるのは、$x _ {t - 1} \rightarrow x _ t$ という遷移が起こる確率密度関数（これは例えば現象に対する微分方程式、差分方程式が利用される事が多い。自己位置推定の文脈ではオドメトリと呼ばれる部分になる）と、$x _ t$ を観測すると $y _ t$ が得られるという観測システムの確率密度関数（これは単に計測機器の誤差を見積もったり、電流値と電圧の状態から電力だけが観測されるなどの物理的変換が記述される）が重要になってきます。
通常、状態空間モデルによって状態遷移モデル $p(x _ t \mid x _ {t-1})$ と観測モデル $p(y _ t \mid x _ t)$ が記述され、その際の観測値による内部状態の推定に逐次ベイズフィルタが利用されるという文脈で利用されます。
余談ですが、俺に加え平滑化（スムージング）という概念も上記の話の流れで導入することができます（ただ、スムージングを利用する例は、データが一通り揃ってから、全体を眺め直すという行為に相当しており、逐次的な処理というわけではないので割愛しました。あと実用上、時系列データが一通り揃ったあとに過去を振り返るような分析は、一般的な統計分析そのものであり、手段としてスムージングはそれほど大きな選択肢に上がってくるものでもない気がしています。無論目的によるのですが…。）
