# Multiple change point model by dynamic programming

**URL:** <https://discourse.mc-stan.org/t/multiple-change-point-model-by-dynamic-programming/1172>\
**Category:** Modeling\
**Created:** [July 7, 2017, 2:23am UTC](https://discourse.mc-stan.org/t/multiple-change-point-model-by-dynamic-programming/1172 "2017-07-07T02:23:33Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [July 7, 2017, 2:23am UTC](https://discourse.mc-stan.org/t/multiple-change-point-model-by-dynamic-programming/1172/1 "2017-07-07T02:23:33Z")

</div>

I am trying to fit a model to time series data on animal step length that allows for three distinct states: early, mid and late. The naming of these states is crucial since in this model the states are temporally-ordered: an animal must go from early-\>mid-\>late (or alternatively remain in early, or transition incompletely from early-\>mid in the course of the experiment); an animal cannot go from early-\>mid-\>early etc.

I use the forward algorithm (see below; I am fitting Speed data that is captured from GPS collars fit to elephants) to estimate the log probability of each possible end state, which I then marginalise and use to increment the total log probability.

My question is – does anyone know how to do the Viterbi algorithm (or its analogue) for the case I describe where the states are ordered? I have had a Google but I cannot seem to find any decent references.

```
data{
  int T;
  int K;
  real Speed[T];
}

transformed data {
  real log_unif;
  int c[2];
  c[1] = 1;
  c[2] = 3;
  log_unif = -log(T);
}

parameters{
  positive_ordered[K] alpha;
}

transformed parameters{
  matrix[T,K] lp;
  
  lp[1] = rep_row_vector(-log(K),K);
  
  for(t in 2:T){
    for(k in 1:K){
      if(k==1) # early (long step length)
        lp[t,k] = lp[t-1,k] + exponential_lpdf(Speed[t]|alpha[k]);
      else if(k==2) # late (medium step length)
        lp[t,k] = log_sum_exp(lp[t-1,2:3]) + exponential_lpdf(Speed[t]|alpha[k]);
      else # mid (short step length during birth of calfs)
        lp[t,k] = log_sum_exp(lp[t-1,c]) + exponential_lpdf(Speed[t]|alpha[k]);
    }
  }
}

model{
  target += log_sum_exp(lp[T]);
  alpha ~ exponential(1);
}
```

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [July 9, 2017, 2:43pm UTC](https://discourse.mc-stan.org/t/multiple-change-point-model-by-dynamic-programming/1172/3 "2017-07-09T14:43:13Z")

</div>

In case this is useful to someone else in the future I have now figured out how to implement Viterbi for the constrained HMM that I mentioned before (see below).

Here because the middle state has a lower step length than the outer two, I reorder the alpha parameter. In my model the state sequence must be 1-\>2-\>3. It cannot go backwards, for example 2-\>1. It also cannot jump from 1-\>3. So from each state the following transitions are permitted:

- State 1: 1-\>1 or 1-\>2.
- State 2: 2-\>2 or 2-\>3.
- Sate 3: 3-\>3 only.

Unlike most HMMs I do not have a transition probability matrix, because really I am just trying to implement a two change point model in an efficient manner.

```
data{
      int T;
      int K;
      real Speed[T];
    }

transformed data {
  real log_unif;
  log_unif = -log(T);
}

parameters{
  positive_ordered[K] alpha;
}

transformed parameters{
  matrix[T,K] lp;
  vector[K] alphaReordered;
  
  alphaReordered[1] = alpha[1];
  alphaReordered[2] = alpha[3];
  alphaReordered[3] = alpha[2];
  
  lp[1] = rep_row_vector(-log(K),K);
  
  for(t in 2:T){
    for(k in 1:K){
      if(k==1) # early (long)
        lp[t,k] = lp[t-1,k] + exponential_lpdf(Speed[t]|alphaReordered[k]);
      else if(k==2) # mid (short)
        lp[t,k] = log_sum_exp(lp[t-1,1:2]) + exponential_lpdf(Speed[t]|alphaReordered[k]);
      else # late (medium)
        lp[t,k] = log_sum_exp(lp[t-1,2:3]) + exponential_lpdf(Speed[t]|alphaReordered[k]);
    }
  }
}

model{
  target += log_sum_exp(lp[T]);
  alpha ~ exponential(1);
}

generated quantities{
  int path[T];
  vector[K] stepLength;
  for(k in 1:K)
    stepLength[k] = 1/alphaReordered[k];
  {
  matrix[T,K] mLogProb;
  int mBestIndex[T-1,K];
  real aLogProbTemp;
  int aBestFinal;
  
  for(t in 1:T){
    for(k in 1:K){
      if(t==1){
        if(k==1)
          mLogProb[t,k] = exponential_lpdf(Speed[t]|alphaReordered[k]);
        else
          mLogProb[t,k] = negative_infinity();
      }else if(t==2){
        if(k==1){
          mLogProb[t,k] = mLogProb[t-1,1] + exponential_lpdf(Speed[t]|alphaReordered[k]);
          mBestIndex[t-1,k] = 1;
        }
        else if(k==2){
          mLogProb[t,k] = mLogProb[t-1,1] + exponential_lpdf(Speed[t]|alphaReordered[k]);
          mBestIndex[t-1,k] = 1;
        }
        else{
          mLogProb[t,k] = negative_infinity();
          mBestIndex[t-1,k] = 0;
        }
      }else if(t==3){
        if(k==1){
          mLogProb[t,k] = mLogProb[t-1,1] + exponential_lpdf(Speed[t]|alphaReordered[k]);
          mBestIndex[t-1,k] = 1;
        }
        else if(k==2){
          mLogProb[t,k] = max(mLogProb[t-1,1:2]) + exponential_lpdf(Speed[t]|alphaReordered[k]);
          aLogProbTemp = negative_infinity();
          for(i in 1:2){
            if(mLogProb[t-1,i]>aLogProbTemp){
              aLogProbTemp = mLogProb[t-1,i];
              mBestIndex[t-1,k] = i;
            }
          }
        }
        else{
          mLogProb[t,k] = mLogProb[t-1,2] + exponential_lpdf(Speed[t]|alphaReordered[k]);
          mBestIndex[t-1,k] = 2;
        }
      }else{
        if(k==1){
          mLogProb[t,k] = mLogProb[t-1,1] + exponential_lpdf(Speed[t]|alphaReordered[k]);
          mBestIndex[t-1,k] = 1;
        }
        else if(k==2){
          mLogProb[t,k] = max(mLogProb[t-1,1:2]) + exponential_lpdf(Speed[t]|alphaReordered[k]);
          aLogProbTemp = negative_infinity();
          for(i in 1:2){
            if(mLogProb[t-1,i]>aLogProbTemp){
              aLogProbTemp = mLogProb[t-1,i];
              mBestIndex[t-1,k] = i;
            }
          }
        }
        else{
          mLogProb[t,k] = max(mLogProb[t-1,2:3]) + exponential_lpdf(Speed[t]|alphaReordered[k]);
          aLogProbTemp = negative_infinity();
          for(i in 2:3){
            if(mLogProb[t-1,i]>aLogProbTemp){
              aLogProbTemp = mLogProb[t-1,i];
              mBestIndex[t-1,k] = i;
            }
          }
        }
      }
    }
  }
  
  aLogProbTemp = negative_infinity();
  for(k in 1:K){
    if(mLogProb[T,k]>aLogProbTemp){
      aLogProbTemp = mLogProb[T,k];
      aBestFinal = k;
    }
  }
  
  for(t in 1:T){
    int aIndex;
    aIndex = T - t + 1;
    if(aIndex == T){
      path[aIndex] = aBestFinal;
    }else if(aIndex==1){
      path[aIndex] = 1;
    }else{
      path[aIndex] = mBestIndex[aIndex,path[aIndex+1]];
    }
  }
  }
}
```

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 13, 2017, 6:57pm UTC](https://discourse.mc-stan.org/t/multiple-change-point-model-by-dynamic-programming/1172/4 "2017-07-13T18:57:38Z")

</div>

Rather than

```
aLogProbTemp = negative_infinity();
for (i in 1:2){
  if (mLogProb[t - 1, i] > aLogProbTemp) {
    aLogProbTemp = mLogProb[t-1, i];
    mBestIndex[t-1, k] = i;
  }
}

```

may I recommend:

```
mBestIndex[t - 1, k] = (mLogProb[t - 1, 1] > mLogProb[t - 1, 2]) ? 1 : 2;

```

It cuts out a dummy (-infinity), looping, and reduces the number of comparisons from one to two.