StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

累積和を使って計算の無駄を省く(変化点検出の例)

メーリングリストでStanにおいて累積和を使って変化点検出を高速化する話がありましたのでメモです。

ここではRにはじめから用意されているNileのデータに対して変化点検出します。プロットすると以下です。

f:id:StatModeling:20201106201449p:plain

ここでは、ある変化点より左の部分では平均mu_l標準偏差sigma正規分布に従い、右の部分では平均mu_r標準偏差sigma正規分布に従うとします。

すると、変化点は離散値をとるパラメータなので、周辺化消去しなくてはいけません。単純にはif_else関数を使った以下の実装になります。

  • 7, 8行目:範囲をおおまかに指定しています。これは実行時にデータを1000で割ってスケーリングするので、この値になっています。

しかし、この実装は各cpにおいて、normal_log(Y[t], mu_l, sigma)normal_log(Y[t], mu_r, sigma)を重複して評価しているので、計算の無駄があります(以下の図)。また、Stanのif_else(cond, x, y)condによらずxyの両方が評価されるのでイマイチです。

f:id:StatModeling:20201106201458p:plain

そこで、累積和を使って実装した例は以下になります。

  • 14,24行目:transformed parametersブロックにおいても途中で{ }でくくっておくと、lp_llp_rはサンプリングされなくなります。tipsです。
  • 23行目:lp_lが黒丸の累積和を保持するvectorに、(lp_r[T] - lp_r)が白丸の累積和を保持するvectorになります(以下の図)。これらを算出するために、21行目と22行目でcumulative_sum関数を使っています。この変更により二重ループが一重ループとなり、計算のオーダーが   \mathcal{O} (T^2) から  \mathcal O (T) に減ります。

f:id:StatModeling:20201106201502p:plain

実行するRコードは以下になります。

  • 3行目:1000で割ってスケーリングしています。
  • 9~10行目:これでどこに変化点があるかの確率分布を算出しています。Stan 2.9.0のマニュアルの「12. Latent Discrete Parameters」の章を参照してください。

計算の無駄を省くことにより、この例題では実行速度はおよそ10倍(僕の環境では9秒→0.9秒)になりました。時点の数が100でこの差が出るので、時点の数がもっと増えたら差は劇的になるでしょう。

ちなみに変化点の確率分布は以下になります。

f:id:StatModeling:20201106201453p:plain

この手の「計算の無駄を省く」方法は競技プログラミング勢は無敵と思いますので、あやかってsucroseさんの記事から通称蟻本のDRMなしのpdfを買ってみました。勉強がてら電車で考えながら読もうと思います。