はじめに
時系列データを扱っていると、株価が下落から上昇に転じた点や、サイトのアクセスカウントが急に増えた点等々のデータの傾向が大きく変わった(レジームシフトした)点を知りたいことがたまにあります。このような"流れ変わったな"という点を見つけるための方法、変化点検知の手法は様々に存在します。
Pythonであれば ruptures や river のようなライブラリがこの分野で特に有名であり、大抵の場合にはこれらのライブラリのいずれかのアルゴリズムを使うことでこの用途を満たす事ができると思います。
しかしその一方、オンラインベイズ変化点検知(Bayesian Online Change Point Detection; BOCPD)は便利な(かつ興味深い)割に実装例が少なかったため、勉強と拡張してみたい欲を兼ねて実装(かなりClaudeのお世話になっていますが
)してみることにしました。
オンラインベイズ変化点検知とは
さて、BOCPDのコアとなるアイデアは「任意の時点で、その直前までどれくらい“同じレジーム(同じ分布から生成されているか)”が続いていたかを確率的に推定する」という点にあります。
この「何回同じレジームが続いている？」というのをBOCPDではランレングス(run length) $r_t$という値の分布で表現します。同じレジームが継続していると推定される場合は尤もらしい $r_t$ の値は増えていき、逆に生成分布が変わりレジームシフトが起きたと推定される場合は $r_t = 0$ の確率値が最も高くなります。
この図はBOCPDの元論文から引用したもので、ランレングスの推移を表すものとなっています。レジームが $g_1$ → $g_2$ → $g_3$ と推移するのに伴ってランレングス $r_t$ がリセットされていることが分かりますね。
多くの変化点検知アルゴリズムでは、レジームシフトが起きたか起きなかったのかのTrue/Falseのみを判定するためのものとなっていますが、一方でBOCPDではレジームシフトが起きたかどうかをランレングスの確率値として継続的に保持しています。そのおかげで、
- レジームシフトが起きていた場合と起きていなかった場合を(その後データを入力していっても)同時に想定できる
- データ点1点ごとに、今まさにレジームシフトが起きているのかを判定できる
といったメリットがあり、そこがBOCPDの興味深い点となっています。
ランレングス分布の逐次更新
では、そのランレングス分布はどのように更新すれば良いのでしょうか？BOCPDでは、新たなデータが入力されるたびに
- 新たなデータが生成分布にとってどのくらいの珍しさか(予測確率)を計算する
- 予測確率を使って、レジームが変化していない確率(成長確率)を計算する
- 予測確率を使って、レジームが変化している確率(変化点確率)を計算する
- 成長確率と変化点確率で規格化して確率値とする
- 新たに観測されたデータを下に生成分布を更新する
のように処理します。「レジームシフトが起きたか？」という条件分岐が全く存在していない点が面白いところだと思います。レジームシフトが起きたかどうかはあくまで確率値として計算されています。
実際に論文中のアルゴリズムを見てみましょう。このアルゴリズムリストの3,4,5が先程の1,2,3に対応しています。
今回実装したプログラムだと以下の部分に対応します(確率値の掛け算をそのまま計算しないよう、対数を取っているので若干異なってます):
数式を見てみると、まず成長確率は
P(r_t = r_{t-1}+1,\; x_{1:t})
=
\underbrace{P(r_{t-1}, x_{1:t-1})}_{\text{$t-1$のランレングス}}
\;
\underbrace{\pi_t^{(r)}}_{\text{新データ$x_t$の予測確率}}
\;
\underbrace{\bigl(1 - H(r_{t-1})\bigr)}_{\text{ランが伸びる想定確率(事前に設定)}}
次に、変化点確率は
P(r_t = 0,\; x_{1:t})
=
\sum_{r_{t-1}}
\underbrace{P(r_{t-1}, x_{1:t-1})}_{\text{$t-1$のランレングス}}
\;
\underbrace{\pi_t^{(r)}}_{\text{新データ$x_t$の予測確率}}
\;
\underbrace{H(r_{t-1})}_{\text{変化点が起こる想定確率(事前に設定)}}
のように計算されています。変化点確率の方のみ$\sum_{r_{t-1}}$で周辺化されているのは、どのようなランレングスの場合からでも変化点が起こり得るから、と考えると直感的ですね。
また、関数 $H(r_{t-1})$ はハザード関数と呼ばれるものであり、この関数によって $r_{t}$ のリセットがどのくらい起こりやすいかを調整しています。 $r_t$ の関数ですので、ランレングスが長くなったほど変化点が起こりやすいという設定もできる一方、一般的には $r_t$ によらず常に一定値をとる設定(無記憶にする)が多いようです。
予測分布について
BOCPDではレジームが継続しているかどうかを生成分布からの予測確率から判断する、というのを前チャプターで(アルゴリズムの3. Evaluate Predictive Probabilityのところです)確認しました。となると、次に気になるのが、その「どのくらい珍しいデータなの？」を計算する、"生成分布"をどうアップデートするのか、という部分かと思います。
BOCPDでは共役事前分布を使うことを想定しています。というのも、新しいデータのたびにMCMCや変分ベイズを行っていると計算コスト的に予測が間に合わなくなってしまう可能性があるためです。
ですので、基本的には事前分布と尤度それぞれに共役事前分布となるペアを設定し、あとはそこから充分統計量としてデータから分布を計算します。それ故にデータの分布に合わせた柔軟なモデリングが可能となっています。
ただし、実際にどのライブラリがどの分布に対応しているのかはまちまちであり、だからこそ今回自分で拡張できるよう、ライブラリを作成したくなったというのが今回の動機です。
実装を見てみる
こちらにコードを置いています。
実装はシンプルさを優先しており、後々に色んな分布を組み込んで遊べるようにしています。ひとまずのところはNormal–Inverse–Gamma事前分布 × ガウス尤度(→t分布予測)の組み合わせを実装しています。
別モデルを追加する際にはgaussian.pyのように抽象クラスの継承からモデルの振る舞いを実装していけばOKです。
作成したモデルはBOCPDクラスのコンストラクタに与えることで使うことができます。
データを入力してみる
では、実際に適当な時系列データを使って変化点検知させてみます。今回はccxtというライブラリを使用してビットコイン相場のチャートを取得し、ボラティリティの変わったところを探してみます。
価格変動のチャートそのものは定常ではないため一日ごとの変動をとり、加えて変化率に注目したいため対数を取った後、BOCPDのモデルにかけてみます。20%以上の確率でレジームシフトが起きたと推定される点をプロットしてみると以下の図のようになりました。
図を見ると、横ばいから上昇/下降へ移行した点、また上昇/下降から横ばいへ移行していそうな点のそれぞれでなんとなく変化点確率が大きくなっていそうなことが分かります。実用するためにはまだまだハザード関数やランレングスの長さの調整が必要そうですが、とりあえずの精度は出ていそうですね！
おわりに
今回はオンラインベイズ変化点検知の考え方を説明し、その実装を行いました。
このBOCPDというアルゴリズムは、事前に想定した分布を新たに観測したデータで更新して予測分布を計算するという、ベイズの中心にあるシステムが自然と組み込まれており、そこに美しさと面白さがあると思います。この記事を通じて少しでもBOCPDの魅力が伝われば幸いです！
おまけ: ほかのベイズ変化点検知ができるPythonライブラリについて
hildensia/bayesian_changepoint_detection がよく使われているようです。
実装だけでなく、ノートブックの使用例が付いているので分かりやすいです！このライブラリでは
| 事前分布 | 尤度 | 予測分布 |
|---|---|---|
| Normal-Gamma分布 | ガウス分布 | Student-t分布 |
| Gaussian-Wishart分布 | 多変量ガウス分布 | 多変量Student-t分布 |
の2通りに対応しています。
また、Facebookが作成しているKatsにもBOCPDが含まれているようです。
ただし、最近はあまりアップデートされていないようでして、私が試したときは「Python3.9じゃないと動かないんだけど？」というissueが立っていました。実際その時にはインストールでコケたような記憶があります。その後のコミットで3.12に対応しているようですね。