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

メーリングリストで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によらずxとyの両方が評価されるのでイマイチです。

f:id:StatModeling:20201106201458p:plain

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

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