Multiple change point model by dynamic programming

Modeling
1

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);
}
3

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]];
    }
  }
  }
}
1 Like
4

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.