Stats Works
  • About This Website

HMMs: Viterbi Algorithm

This investigates HMMs (Hidden Markov Models) with two possible states and six possible emissions using the viterbi algorithm.

Given a Markov process with two states (state 1 and state 2) and an arbitrary number of possible observed emissions , with transition probabilities and emission probabilities given for each state, we can find the most likely sequence of states given some observed emissions using the viterbi algorithm. This algorithm is an example of dynamic programming which produces results far faster than a brute force approach.

This takes less time to run than considering all possible sequences of states -- we only need to consider the state probabilities forwards, rejecting paths that are unlikely at each point.

The example which the following function defaults to is the case of a fair and loaded die, each with some probability of values 1-6. The roller does not show you what die he is using, but we know that he will generally switch between the two die with some transition probs. The goal is then to find the most likely sequence of fair and loaded die given the observed sequence of emissions (for example 1,2,2,1,3).

In [12]:
viterbialgorithm <- function(values, initprobstates = c(.5,.5), 
                             probvectorstate1 = c(rep(1/6,6)), probvectorstate2 = c(rep(1/10,5),1/2), 
                             Ps1s2 = 0.05, Ps2s1 = .1){
  ## values (vector): a series of observed values of two emissions.
  ## initprobstates (vector, len 2): the probabilities of entering the system in either state. (assuming .5, .5 is the least informative)
  ## probvectorstate1 (vector): the emission probabilities of the first state.
  ## probvectorstate2 (vector): the emission probabilities of the second state. 
  ## Ps1s2 (numeric): probability of entering state 2 from state 1. 
  ## Ps2s1 (numeric): probability of entering state 1 from state 2. 
  
  
  iterations <- length(values)
  
  initprobstate1 = initprobstates[1]
  initprobstate2 = initprobstates[2]
  
  P11 <- 1 - Ps1s2
  P22 <- 1 - Ps2s1
  
  if (iterations < 1) {
    stop("Length of vector must be a positive integer.")
  }
  
  # Initializes the ret list for the first sequence. 
  ret <- list(c(s0 = 0,s1 = -Inf, s2 = -Inf, st = 0, prevst = NA))
  
  ## Creates the transition matrix. 
  transmatrix <- matrix(c(0,0,0,initprobstate1,P11,Ps2s1,initprobstate2,Ps1s2,P22),ncol = 3)
  
  maxfun <- function(potentialstate, prevvector){
    ## potentialstate (integer): integer referring to the potential state (1 or 2)
    ## prevvector (vector): previous probability of states. 
    
    trans <- transmatrix[, potentialstate + 1]
    mx <- max(prevvector[1:3] + log(trans))
    prevst <- which.max(prevvector[1:3] + log(trans)) - 1
    return(list(val = mx,state = prevst))
  }
  
  for (i in 1:iterations) {
    maxst1 <- maxfun(1,ret[[i]])
    maxst2 <- maxfun(2,ret[[i]])
    ekst1 <- log(probvectorstate1[values[i]])
    ekst2 <- log(probvectorstate2[values[i]])
    
    V1 <- maxst1[['val']] + ekst1
    V2 <- maxst2[['val']] + ekst2
    
    
    end <- which.max(c(V1,V2)) 
    ret[[i + 1]] <- c(s0 = -Inf, s1 = V1, s2 = V2, st = end, prevst = c(maxst1[['state']],maxst2[['state']]))
  }
  
  finalstate <- ret[[iterations + 1]]['st']
  path <- c(rep(0,iterations + 1))
  path[iterations + 1] <- finalstate
  
  for (i in (iterations):2) {
    path[i] <- ret[[i + 1]][5:6][path[i + 1]]
  }
  
  res <- as.data.frame(ret[-1])
  colnames(res) <- c(values)
  
  return(list(steps = res,path = path))
}
In [13]:
observed <- c(6,5,1,1,6,6,4,5,
              3,1,3,2,6,5,1,2,
              4,5,3,6,6,6,4,6,
              3,1,6,3,6,6,6,3,
              1,6,2,3,2,6,4,5,
              5,2,3,6,2,6,6,6,
              6,6,6,2,5,1,5,1,
              6,3,1)

results <- viterbialgorithm(observed)
resultswords <- c("Fair","Loaded")[results$path]
resultswords
  1. 'Fair'
  2. 'Fair'
  3. 'Fair'
  4. 'Fair'
  5. 'Fair'
  6. 'Fair'
  7. 'Fair'
  8. 'Fair'
  9. 'Fair'
  10. 'Fair'
  11. 'Fair'
  12. 'Fair'
  13. 'Fair'
  14. 'Fair'
  15. 'Fair'
  16. 'Fair'
  17. 'Fair'
  18. 'Fair'
  19. 'Fair'
  20. 'Loaded'
  21. 'Loaded'
  22. 'Loaded'
  23. 'Loaded'
  24. 'Loaded'
  25. 'Loaded'
  26. 'Loaded'
  27. 'Loaded'
  28. 'Loaded'
  29. 'Loaded'
  30. 'Loaded'
  31. 'Loaded'
  32. 'Loaded'
  33. 'Loaded'
  34. 'Loaded'
  35. 'Loaded'
  36. 'Loaded'
  37. 'Loaded'
  38. 'Loaded'
  39. 'Loaded'
  40. 'Loaded'
  41. 'Loaded'
  42. 'Loaded'
  43. 'Loaded'
  44. 'Loaded'
  45. 'Loaded'
  46. 'Loaded'
  47. 'Loaded'
  48. 'Loaded'
  49. 'Loaded'
  50. 'Loaded'
  51. 'Loaded'
  52. 'Fair'
  53. 'Fair'
  54. 'Fair'
  55. 'Fair'
  56. 'Fair'
  57. 'Fair'
  58. 'Fair'
  59. 'Fair'

So this tells us that given the observed die rolls, the most likely sequence of states is as shown above.

This did not have to be done manually -- we can use the HMM package to do this. It will serve as a check.

In [14]:
require(HMM)
mod <- initHMM(c("Fair","Loaded"),c(1,2,3,4,5,6),startProbs = c(.5,.5),transProbs = rbind(c(.95,.05),c(.1,.90)),emissionProbs = rbind(c(rep(1/6,6)),c(rep(1/10,5),1/2)))
viterbi(mod,observed) == resultswords
## :)
  1. TRUE
  2. TRUE
  3. TRUE
  4. TRUE
  5. TRUE
  6. TRUE
  7. TRUE
  8. TRUE
  9. TRUE
  10. TRUE
  11. TRUE
  12. TRUE
  13. TRUE
  14. TRUE
  15. TRUE
  16. TRUE
  17. TRUE
  18. TRUE
  19. TRUE
  20. TRUE
  21. TRUE
  22. TRUE
  23. TRUE
  24. TRUE
  25. TRUE
  26. TRUE
  27. TRUE
  28. TRUE
  29. TRUE
  30. TRUE
  31. TRUE
  32. TRUE
  33. TRUE
  34. TRUE
  35. TRUE
  36. TRUE
  37. TRUE
  38. TRUE
  39. TRUE
  40. TRUE
  41. TRUE
  42. TRUE
  43. TRUE
  44. TRUE
  45. TRUE
  46. TRUE
  47. TRUE
  48. TRUE
  49. TRUE
  50. TRUE
  51. TRUE
  52. TRUE
  53. TRUE
  54. TRUE
  55. TRUE
  56. TRUE
  57. TRUE
  58. TRUE
  59. TRUE
Sean Ammirati - creator of Stats Works. He can be reached on Github, LinkedIn and email.
Comments
comments powered by Disqus

Articles in Markov Chains

  • Hidden Markov Models Extensions »

All Articles

  • « Hotelling's T-Test Example
  • Constrained Least Squares »

Published

Dec 9, 2018

Category

Markov Chains

Tags

  • HMMS 2
  • markov chains 2
  • viterbi 1

Contact

  • Stats Works - Making Statistics Work for You