# Identify change points in multiple time-series

**URL:** <https://discourse.mc-stan.org/t/identify-change-points-in-multiple-time-series/8999>\
**Category:** Modeling\
**Tags:** techniques\
**Created:** [May 26, 2019, 7:47pm UTC](https://discourse.mc-stan.org/t/identify-change-points-in-multiple-time-series/8999 "2019-05-26T19:47:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![fabio](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fabio/32/935_2.png) [@fabio](https://discourse.mc-stan.org/u/fabio)\
**Post date:** [May 26, 2019, 7:47pm UTC](https://discourse.mc-stan.org/t/identify-change-points-in-multiple-time-series/8999/1 "2019-05-26T19:47:14Z")

</div>

I have l = 1 \dots L, arbitarily long time-series, each one with Nt\_l normally distributed measures D\_l .

Each serie contains a change-point located at an arbitrary time \tau\_l , so that:

D\_l(t) \sim \begin{cases} \mathcal{N}(\mu1\_{p,l}, \sigma\_p) & \mbox{if } t \lt \tau\_l \\ \mathcal{N}(\mu2\_{p,l}, \sigma\_p), & \mbox{if } t \geq \tau\_l \end{cases}

Moreover, it is known that \mu1\_{p,l} \sim \mathcal{N}(mu1, \sigma) and \mu2\_{p,l} \sim \mathcal{N}(mu2, \sigma).

Since I would write a Stan model to do inference on my measures, firstly I generated some synthetic data to check if the parameters (especially the \tau) are correctly calculated:

```
L = 3
Nt = 50

tau = sample(x = seq(from = round(Nt /3), to = round(2 * Nt / 3)), size = L, replace = TRUE)
toEnd = Nt - tau

mu1 = -0.7
mu2 = 1.3
sigma = 0.4

mu1_p = rnorm(L, mean = mu1, sd = sigma)
mu2_p = rnorm(L, mean = mu2, sd = sigma)

sigma_p = 0.2

D = matrix(data = NA_real_, nrow = Nt, ncol = L)

for (t in 1: Nt) {
    for (l in 1: L) {
        D[t, l] = ifelse(t < tau[l], rnorm(1, mean = mu1_p[l], sd = sigma_p), rnorm(1, mean = mu2_p[l], sd = sigma_p))
    }
}

```

In this case, all the L=3 time-series have the _same size_ Nt=50 which it is **not** the real case I will apply the model to. This comes to the fact that I started coding the Stan model extending the [nowave.it](http://nowave.it/pages/bayesian-changepoint-detection-with-r-and-stan.html)’ script adapting it to the model above. The consequence is that D id a 50x3 matrix. The three \tau are generated to lay somewhat in the middle of each time-series, and this is consistent with the physical mechanism I want to describe.  
So, here’s the translation of the math model to Stan:

```stan
data {
    int<lower=1> Nt;
    int<lower=1> L;
    vector[L] D[Nt]; 
}

parameters {
    real mu1;
    real mu2;
    real<lower = 0> sigma;
    vector[L] mu1_p_raw;
    vector[L] mu2_p_raw;
    real<lower = 0> sigma_p;
}

transformed parameters {
    vector[L] mu1_p = mu1 + sigma * mu1_p_raw;
    vector[L] mu2_p = mu2 + sigma * mu2_p_raw;
    vector[Nt] log_p[L];
    
    for (l in 1:L) {
        vector[Nt + 1] log_p_e;
        vector[Nt + 1] log_p_l;
        log_p_e[1] = 0;
        log_p_l[1] = 0;
        for (i in 1:Nt) {
            log_p_e[i + 1] = log_p_e[i] + normal_lpdf(D[i, l] | mu1_p[l], sigma_p);
            log_p_l[i + 1] = log_p_l[i] + normal_lpdf(D[i, l] | mu2_p[l], sigma_p);
        }
        log_p[l] = rep_vector(-log(Nt) + log_p_l[Nt + 1], Nt) + head(log_p_e, Nt) - head(log_p_l, Nt);
    }
}

model {
    mu1 ~ std_normal();
    mu2 ~ std_normal();
    
    sigma ~ std_normal();
    
    mu1_p_raw ~ std_normal();
    mu2_p_raw ~ std_normal();
    
    sigma_p ~ std_normal();
    
    for (l in 1:L)
        target += log_sum_exp(log_p[l]);
}

generated quantities {
    int<lower = 1, upper = Nt> tau[L];
    for (l in 1:L) {
        tau[l] = categorical_rng(softmax(log_p[l]));
    }
    
}

```

I am quite happy with the results, although the reparametrization of \mu1\_p and \mu2\_p and `adapt_delta=0.99` do not have eliminated some sporadic divergent transitions after warmup. Here are the results:

```
               mean se_mean sd 50% n_eff Rhat
mu1_p[1] -1.3634953 0.0005508518 0.03529051 -1.3634321 4104.3675 0.9999053
mu1_p[2] -0.2586235 0.0005489039 0.03590872 -0.2591301 4279.6400 0.9993275
mu1_p[3] -0.7160250 0.0007111709 0.04712015 -0.7162227 4390.0093 1.0004151
mu2_p[1] 1.2859460 0.0006491305 0.04319718 1.2852867 4428.3997 0.9994939
mu2_p[2] 1.3433570 0.0007406048 0.04210413 1.3435833 3232.0374 1.0004466
mu2_p[3] 1.8209715 0.0005387842 0.03373185 1.8207860 3919.6797 0.9991053
sigma_p 0.1962541 0.0002211053 0.01155979 0.1961298 2733.3879 1.0001940
mu1 -0.7018092 0.0092852402 0.33567857 -0.7159256 1306.9559 1.0012539
mu2 1.3134600 0.0112513183 0.37241814 1.3493183 1095.6071 1.0013171
sigma 0.5961289 0.0080242368 0.25166286 0.5370761 983.6277 1.0025183
tau[1] 31.0000000 NaN 0.00000000 31.0000000 NaN NaN
tau[2] 30.0000000 NaN 0.00000000 30.0000000 NaN NaN
tau[3] 19.0000000 NaN 0.00000000 19.0000000 NaN NaN

```

wich are consistent with the original parameters:

```
> mu1_p
[1] -1.3758223 -0.2042016 -0.7435864
> mu2_p
[1] 1.253103 1.373233 1.812222
> tau
[1] 31 30 19

```

So far, so good, but this is **not** the structure of my measures because I should have arbitrarily long time series. The input structure should be similar to:

```
L = 3
Nt = sample(x= seq(from = 50, to = 100), size = L, replace = TRUE)

tau = c()
for (l in 1:L) {
    tau = c(tau, sample(x = seq(from = round(Nt[l] /3), to = round(2 * Nt[l] / 3)), size = 1))
}

N = sum(Nt)

D = data.frame()
for (l in 1: L) {
    for (n in 1: Nt[l]) {
        D = rbind(D, data.frame(time= n, IDx = l, D = ifelse(n < tau[l], rnorm(1, mean = mu1_p[l], sd = sigma_p), rnorm(1, mean = mu2_p[l], sd = sigma_p))))
    }
}

```

`D` now becomes a `data frame` and each series is identified by the value of `IDx`.

And here comes the question:  
With this structure I am not able (in Stan) to correctly indexing the `vector[N] log_p;` in the `model` block (`log_p` is no more a matrix, since the series are not equally long) and indexing when marginalizing out the \tau in the `generated quantities` block.

Any help is appreciated.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [May 30, 2019, 9:24am UTC](https://discourse.mc-stan.org/t/identify-change-points-in-multiple-time-series/8999/2 "2019-05-30T09:24:40Z")

</div>

Since nobody replied, I will give it a try.

But first:

> [@fabio](#):
>
> although the reparametrization of μ1p\mu1\_p and μ2p\mu2\_p and `adapt_delta=0.99` do not have eliminated some sporadic divergent transitions after warmup

It usually does not pay off to ignore those as they often bite you later. I don’t deeply understand changepoint models, so can’t help you directly, but I’ve written about a few tricks to handle divergences at [Taming Divergences in Stan Models](https://www.martinmodrak.cz/2018/02/19/taming-divergences-in-stan-models/)

But on to the implementation question:  
What I would do would be to concatenate all the time series into a single vector and manually unserialize them as needed. It could look something like (not tested, just a sketch, but I made it work for a different task):

```stan
data {
  int<lower=1> Nt;
  int<lower=1> lengths;
  vector[sum(lengths)] D;
}

transformed data {
  int<lower=0> series_start[Nt];

    series_start[1] = 1;
    for(i in 2:Nt) {
      series_start[i] = series_start[i - 1] + lengths[i - 1];
    }

}

transformed parameters {
  vector[sum(lengths)] log_p;
  for(i in 1:Nt) {
     vector[lengths[i]] log_p_for_series;
     //compute log_p for the given series
     //...
     log_p[series_start[i] : (seris_start[i] + lengths[i])] = log_p_for_series;
  }
}

```

Hope that helps and best of luck!