2007年に発表された、Bayesian Online Changepoint Detectionの紹介と解説を行います。該当論文のarxivへのリンクはこちらになっています。この解説記事は、この論文の英語版の解説記事を参考にしています。
Bayesian Online Changepoint Detection は変化点検地に関する論文であり、時系列データなどのデータが与えられたとき、データ列の傾向が変化するような点を検知するようなアルゴリズムです。時系列データで例えると、例えば株価予測をしていて、今まで上昇トレンドだった株価が急に下落したようなタイミングで変化点を検知するようなアルゴリズムになっています。
手法としてはベイズ推論の枠組みを利用しています。この論文は非常に有名で数多くの引用がされているものなので、一度読んだことがある人も多いのではないでしょうか。今回はこのBayesian Online Changepoint Detectionの紹介と解説をしていきます。
上の図(a)の$X_t$が、時刻$t$における時系列データの値だと考えてください。横軸が$t$であり、時系列データを考えている場合は時刻$t$を示しています。
Bayesian Online Changepoint Detectionは基本的に、様々なデータで扱える汎用的なアルゴリズムですが、今回は時系列データの変化点検知を例に説明をします。ここで、図（a）の時系列を見ると、パッと見て時系列データが3つのパートに分けられるのがわかると思います。
図(a)の中でも、$g_1, g_2, g_3$に分けられており、それぞれのパートの間で点線が引かれているのがわかると思います。さらに下の図を見ると、$r_t$という変数が導入されており、これが前のパートから何データ分経過したかを示す領となっています。
今回のこのBaysian Change Poiint Detectionで予測するのは、この$r_t$であり、$r_t = 0$となった時、パートの変化点を検知したというようになります。（解釈としては、前の変化パートから0つ経過した、つまり、今回変化点が訪れてカウンタがリセットされた、ということになります）
つまり、$r_t$については下記のようになります。
\begin{equation}
r_t=
\begin{cases}
0 & \text{if \space 時刻tがchange point(変化点)だった場合} \\
r_{t-1}+1 & \text{if \space それ以外}
\end{cases}
\end{equation}この$r_t$のことと、元々の論文中ではrun length と表現しています。run length は、change pointが来ていれば0になり、それ以外であれば1つ増えるという変数になります。
今回、このBayesian ChangePoint Detectionで予測したい確率は$p(r_t | r_{t-1})$になります。
ここでは、$r_{t-1}$ が起こった際の、$p(r_t)$ の条件付き確率になります。 さらにここから、一番求めたいのは、時刻1~$t$までの全データ$X_{1:t}$が与えられた時の、$p(r_t | X_{1:t})$と、事後予測分布である、$p(x_{t+1} | X_{1:t})$になります。
ベイズ推定の復習
まず、このアルゴリズムを理解する前に、ベイズ推定を復習しておきましょう。ベイズ推定といえば、あるモデルにおいて、データが与えられた時に、そのモデルのパラメータ$\theta$も確率変数として捉え、定式化するものです。ここら辺は十分にわかっている人は、このベイズ推定の部分は飛ばしても大丈夫です。
ベイズの予測分布とは
さらに、ベイズ推定では、パラメータを確率変数として捉え定式化するだけでなく、確率変数としたパラメタをもとで、次に来るデータの確率分布を導出することができます。これをベイズ推定における予測分布と言います。
これは数式で表現すると、$p(\hat{x} | \theta)$となり、下の数式で表現することができます。ちなみに、この$p(\hat{x} | \theta)$を事後予測分布(posterior predictive distribution)と呼んだりします。
\begin{equation}
p( \hat{x} | \mathcal{D}) = \int p (\hat{x} | \mathcal{D}, \theta) p( \theta |\mathcal{D}) d\theta
\end{equation}
確率変数$x$の取りうる値が確率分布として出てきます。数式では予測分布と言われてイメージできるかもしれませんが、上のようなグラフ（これが確率分布ですが）が得られるとすれば、これは結構嬉しいですよね。
(2)式の右辺の意味は少し難しいです。パラメータ$\theta$で周辺化を行うことで、左辺の予測分布を得ています。このあたりの解説は、かなり長くなるので、別の記事で解説することにします
ベイジアン変化点検出の定式化
さて、ベイズの予測分布の式を踏まえて、今回解説するベイジアン変化点検出の論文の解説に戻ります。今回、所与のデータ$X_{1:t}$がある時の、$x_{t+1}$の予測分布$p(x_{t+1} | X_{1:t})$を考えます。
\begin{equation}
p(x_{t+1} | X_{1:t}) = \sum_{r_t} p(x_{t+1} , r_t | X_{1:t})
\end{equation}まず、考えたい予測分布$p(x_{t+1} | X_{1:t})$ は、(3)のように式変形できます。これは、$p(x_{t+1})$を、パラメータ$r_t$で周辺化(marginalization)しています。
(3)の式変形がピンとこない人のために、かなり簡易化した説明をします。
今完結のため、$X$の取りうる値は、0, 1, 2,3 の4種類で、 $r_t$の取りうる値が、0, 1, 2 の3種類であるという風にすると、下記のようなマトリクスになりますよね。
| $r_t = 0$ | $r_t = 1$ | $r_t = 2$ | 合計($r_t$で周辺化) | |
| $x = 1$ | $\frac{1}{12}$ | $\frac{1}{4}$ | $\frac{1}{6}$ | $\frac{1}{2}$ |
| $x = 2$ | $\frac{1}{18}$ | $\frac{1}{18}$ | $\frac{1}{18}$ | $\frac{1}{6}$ |
| $x = 3$ | $\frac{1}{9}$ | $\frac{1}{9}$ | $\frac{1}{9}$ | $\frac{1}{3}$ |
ここで、(3)式にまた戻ります。さらに(3)式はこのように変形することができます。
\begin{equation}
\begin{split}
p(x_{t+1} | X_{1:t}) &= \sum_{r_t} p(x_{t+1} , r_t | X_{1:t}) \\
&=
\sum_{r_t}
\underbrace{p(x_{t+1} | r_t, \bm{x}_t^{(r)}) }_{\color{red}{UPM \space predictive}}
\space
\underbrace{ p ( r_t | X_{1:t}) }_{\color{blue}{RL \space posterior}}
\end{split}
\end{equation}こうこの辺でかなり、がちゃがちゃしてきました…。
新しい$\bm{x}_t^{(r)}$ は、run length が続いた時の Xの値の集合になります。 例えば、時刻 $t$において、run length が4だった時、$\bm{x} = X_{t-r+1:t}$となります。
ここら辺は、かなりわかりにくですね…。とはいえ、この(4)式がこの論文の手法を理解する上で最初の出発点になるところなので、頑張って理解しましょう。
ここで、冒頭で紹介した Gregorygundersenの解説記事では、(4)式の二行目に登場するそれぞれの確率をそれぞれ UPM (underlying probabilistic model) と、RL posterior (RL事後分布)と表現しています。
ここで、まず、RL 事後分布を考えていきましょう。(UPM 予測分布についても、その次に深掘っていきます）
\begin{equation}
p ( r_t | X_{1:t}) = \frac {p( r_t , X_{1:t}) }{p(X_{1:t})}
\end{equation}が成り立ちます。ここで、(5)式右辺の分子に着目すると、このように式変形することができます。
\begin{equation}
\begin{split}
p( r_t , X_{1:t}) &= \sum_{r_{t-1}} p( r_t, r_{t-1}, X_{1:t}) \\
&= \sum_{r_{t-1}} p( r_t, x_t | r_{t-1}, X_{1:t-1} ) \space p( r_{t-1}, X_{1:t-1}) \\
&= \sum_{r_{t-1}} p(r_{t} | r_{t-1}, X_{1:t-1}) p(x_t | r_t, r_{t-1}, \bm{x}_t^{(r)}) p(r_{t-1} | X_{1:t-1}) \\
&= \sum_{r_{t-1}} p(r_{t} | r_{t-1}) p(x_t | r_{t-1}, \bm{x}_t^{(r)}) p(r_{t-1} | X_{1:t-1})
\end{split}
\end{equation}この式変形は詳しく追っていく必要性があります。元の論文の中でもかなり説明は省いています。
最初の1行目は、確率変数$r_{t-1}$に対して周辺化をしているだけです。次の2行目は、条件付き分布の連鎖律(chain rule)の公式を利用しています。
最後の式変形は、本手法のヒューリスティックな部分が入っているところになるので、丁寧に追っていきます。
p(r_t | r_{t-1}, \bm{x}_{1:t-1}) = p(r_t | r_{t-1}) \\
p(x_t | r_t, r_{t-1}, \bm{x}_{1:t-1}) = p(x_t | r_t, \bm{x}_t^{(r)})このような、ヒューリスティック（仮定）を置いていることに注意しましょう。
UPM predictive の計算
ここで、(4)式における、UPM predictive (UPM予測分布)を深掘っていきましょう。
p(x_{t+1} | r_t, \bm{x}_t^{(l)}) = \int p(x_{t+1 } | \eta ) p(\eta |~ r_t, \bm{x}_{1:t}) d \eta