累積和を使った周辺化消去を用いたベイズ推定による変化点検出の試み（その３）

##1.はじめに
[KFASパッケージを使ったガウス状態空間モデルによる変化点検出の試み](https://qiita.com/kenkenvw/items/68777b8fc5122e9258b0)
[状態空間モデルを使った変化点検出の試み（その２）](https://qiita.com/kenkenvw/items/3fa3c313e372db1d375b)
今回は状態空間モデルではなく、変化点=cp(単位:day)を用いてRstanによる直接的に変化点検出を試みます。
変化点cpは離散値のため、パラメータの推定に離散値が直接扱える、例えばPyMCだと比較的簡単ですが、stanは連続値しか扱えないので工夫が必要になります。

##2.データの準備
ナイル川の流量のデータを準備します。
データは最大値で割ることで正規化しておきます。

```{r}
library(rstan)
library(tidyverse)
library(bayesplot)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

df = data.frame(Time=1:length(Nile), Y=as.vector(Nile)/max(Nile))
data = list(T=length(Nile), Y=as.vector(Nile)/max(Nile))
ggplot(df,aes(x=Time,y=Y)) + geom_line() 
```
![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/6bbc91bb-0e89-e16d-81e7-f451c134b926.png)

##3.stan codeの準備

以下、stan_codeです。
ファイル名は、model.stanで保存します。

```
data {
  int T;
  vector[T] Y;
}

parameters {
  real<lower=0,upper=1> mu_l;
  real<lower=0,upper=1> mu_r;
  real<lower=0> sigma;
}

transformed parameters {
  vector[T] lp;
  vector[T] lp_l;
  vector[T] lp_r;
    for (cp in 1:T) {　　　　
      lp_l[cp] = normal_lpdf(Y[cp] | mu_l, sigma);
      lp_r[cp] = normal_lpdf(Y[cp] | mu_r, sigma);
    }
    lp_l = cumulative_sum(lp_l);
    lp_r = cumulative_sum(lp_r);
    lp = lp_l + (lp_r[T] - lp_r);
  
}

model {
  target += log_sum_exp(lp);
}
```


###3.1対数尤度の計算
平均パラメータは左側と右側の２つ（mu_l,mu_r）を設定しています。
正規化して0から1までしか範囲をとらないため、パラメータも、0〜1に制限しています。


lpの計算は、変化点はcpでmu_lからmu_rに切り替わるため、直接的には次のコード１のとおりです。

```
コード１
 vector[T] lp;
  lp = rep_vector(0, T);
  for (cp in 1:T)
    for (t in 1:T)
      lp[cp] = lp[cp] + normal_lpdf(Y[t] | if_else(t <= cp, mu_l, mu_r), sigma);
```

変化点cpはlpの中に周辺化消去されてなくなります。
しかし、Tについて２重ループするため、計算速度が遅くなります。

コード２は、cumulative_sum（累積和）をとることで、２重ループを解消しています。こうすることで、計算スピードが早くなります。

```
コード２
 for (cp in 1:T) {　　　　
      lp_l[cp] = normal_lpdf(Y[cp] | mu_l, sigma);
      lp_r[cp] = normal_lpdf(Y[cp] | mu_r, sigma);
    }
lp_l = cumulative_sum(lp_l);
lp_r = cumulative_sum(lp_r);
lp = lp_l + (lp_r[T] - lp_r);
```
![logsum.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/48e81f77-9755-8e3e-b795-ea64cd66c30f.png)
lp_lとlp_rの計算方法を図で示しています。

###3.2stanの実行

```{r}
fit = stan(file="model.stan",data=data)

$> 省略
$>SAMPLING FOR MODEL 'model' NOW (CHAIN 4).
$>Chain 4: 
$>Chain 4: Gradient evaluation took 4.3e-05 seconds
$>Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
$>Chain 4: Adjust your expectations accordingly!
$>Chain 4: 
$>Chain 4: 
$>Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
$>Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
$>Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
$>Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
$>Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
$>Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
$>Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
$>Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
$>Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
$>Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
$>Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
$>Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
$>Chain 4: 
$>Chain 4:  Elapsed Time: 0.125205 seconds (Warm-up)
$>Chain 4:                0.131506 seconds (Sampling)
$>Chain 4:                0.256711 seconds (Total)
$>Chain 4: 
```
収束を見てみます。

```{r}
mcmc_rhat(rhat(fit))
mcmc_combo(fit, pars = c("mu_l", "mu_r"))
```
![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/b25c82f7-1c15-5045-efdd-1e8d940e56fe.png)
![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/a3df3224-e70d-c64e-a673-f3e1cff20543.png)
Rhatは1と綺麗に収束しています。
mu_lとmu_rの事後分布のトレースも４chainがきれいに混ざっています。

###4.結果の解釈

```{r}
fit

$>Inference for Stan model: model.
$>4 chains, each with iter=2000; warmup=1000; thin=1; 
$>post-warmup draws per chain=1000, total post-warmup draws=4000.
$>
$>            mean se_mean    sd    2.5%    25%    50%    75% 97.5%
$>mu_l        0.80    0.00  0.02    0.77   0.79   0.80   0.81  0.84
$>mu_r        0.62    0.00  0.01    0.60   0.61   0.62   0.63  0.64
$>sigma       0.10    0.00  0.01    0.08   0.09   0.09   0.10  0.11
$>lp[1]      45.96    0.15  9.49   23.91  40.41  46.99  52.97 61.04
$>lp[2]      48.70    0.14  8.99   27.74  43.47  49.70  55.38 62.90
$>  省略
```
lp[i]はi時点における２つの平均（"mu_l", "mu_r"）の混合対数尤度の事後分布を示します。lpが最も高い時点(i)が変化点としての確率が最も高いところになります。

lpの中央値を推定値として、確率を計算します。

```{r}
fit2 = rstan::extract(fit)
d2 = exp(apply(fit2,2,median))
d3 = data.frame(prob = d2/sum(d2))
d3$time = 1:nrow(d3)
which.max(d3$prob)
ggplot(d3,aes(x = time,y = prob)) + geom_bar(stat="identity")

$>[1] 28
```
![04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/00093b48-0ba4-0b07-cbd9-092ab8b59087.png)
28時点で確率が最も高くなりました。
変化する前の最大時点を推定しているため、次の時点（+1）が変化点になります。
前回の解析の29時点（1899年）と同じになりました。

流量の変化の50%中央値と95%信用区間を計算します。

```{r}
quantile(ms$mu_l,prob = c(0.025,0.5,0.975))*max(Nile)
quantile(ms$mu_r,prob = c(0.025,0.5,0.975))*max(Nile)

$>    2.5%      50%    97.5% 
$>1049.729 1097.622 1146.958 
$>    2.5%      50%    97.5% 
$>820.3091 850.3095 881.5139
```
流量は、変化点において1097.6から850.3に変化しました。

##5.参考

全面的に参考にしました。
StatModeling Memorandum
[累積和を使って計算の無駄を省く（変化点検出の例）](https://statmodeling.hatenablog.com/entry/cumulative_sum-to-reduce-calculation)