Using Hidden Markov Models and Viterbi Algorithm to predict your friend’s moods

HMMs Predict Ramiz’ Mood

• Ramiz’ mood can be either angry or calm
• His mood affects the food he chooses to eat for dinner
• His dinner can be either pizza, burger, or tacos
• If you tell me the food Ramiz chose to eat every day for a month, we can figure out whether Ramiz was angry or calm.
• In HMMs, Ramiz’ mood is a hidden state, while the dinner he eats are the observations.

For simplicity, we will designate food with numbers.

• Pizza = 3
• Burger = 2
• Tacos = 1

We will designate the hidden states (moods) with letters

• Calm = A
• Angry = B

Load the environment

library(tidyverse)
library(HMM)
set.seed(123)

First, lets take a look at the food Ramiz ate last month

data$Visible

##  [1] 1 2 2 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 2 3 2 3 3 2 1 3 2 3

First lets try to generate what we think would be a good Hidden Markov Model for Ramiz’ mood and diet.

We can assume both of his moods are equally likely. This translates as equal priors.

priors = c(
  A=0.5,
  B=0.5
)

We can assume that if he is angry today, he is equally likely to be angry or calm tomorrow. Conversely, if he is calm today, he is equally likely to be angry or calm tomorrow. This would mean he has equal transition probabilities.

# transition probabilities
transition.probabilities = matrix(
  nrow = 2,
  ncol = 2,
  data=rep(0.5, 4)
)
colnames(transition.probabilities) = c('A', 'B')
rownames(transition.probabilities) = c('A', 'B')

Now we have to define the emission probabilities. These are the probabilities that Ramiz will eat a certain type of food for dinner given his mood.

From our time being Ramiz’ friend, we have observed that he is more likely to eat pizza than burgers, and more likely to eat bugers than tacos. We also think he is slightly more likely to eat tacos when he is angry.

# emission probabilities
emission.probabilities = data.frame(
  A = c(
    0.111,
    0.333,
    0.555
  ),
  B = c(
    0.166,
    0.333,
    0.500
  )
)

Finally, let’s save Ramiz’ previous 31 dinner in a vector

# visible/observed data
visible=data$Visible

OK! We have everything we need to define our Hidden Markov Model. But how do we use the model to determine Ramiz’ moods, aka, his hidden states? One popular way to solve this is to use the Viterbi Algorithm. Let’s take a look at an example implementation of this algorithm. It will use as input the priors, emissions and transition probabilities we have defined previously, as well as with the vector of observed states.

# calculate viterbi
My.Viterbi = function(
  emission.probabilities = emission.probabilities,
  transition.probabilities = transition.probabilities,
  priors = priors,
  visible = visible
)
{
  # get number of observations and hidden states
  num.obs=length(visible)
  num.hidden=length(priors)
  
  # get the names for the hidden states
  hidden.states=colnames(emission.probabilities)
  
  # initiate a list object to store the hidden state probabilities
  hidden.states.list=list()
  
  # go through observed states and get largest probabilities
  # for each hidden state
  for (i in 1:num.obs) 
  {
    # get current observation
    curr.visible=visible[i]
    
    # initiate a vetor to store the probabilities
    p.hidden=c()
    
    # get all max probabilities for each hidden state
    # we will work in the logarithm space to avoid
    # numerical problems
    for (h in 1:num.hidden) 
    {
      # get current hidden state
      curr.hidden.state=hidden.states[h]
      
      # get prior prob 
      prior.p=priors[h]
      prior.p=log(prior.p)
      
      # get emission probability
      curr.emission.p=emission.probabilities[curr.visible,h]
      curr.emission.p=log(curr.emission.p)
      
      # get transition probabilities
      if(i>1)
      {
        # probabilities for previous states
        prev.times = hidden.states.list[[i-1]]
        
        # identify the correct transition probabilities
        curr.transition.ps=transition.probabilities[,h]
        curr.transition.ps=log(curr.transition.ps)
        
        # get final transition probabilities
        final.transition.ps=curr.transition.ps+prev.times
      } else 
      {
        final.transition.ps=rep(0, num.hidden)
      }
      
      # get final probabilities for this hidden state
      final.curr.prob=prior.p+final.transition.ps+curr.emission.p
      
      # max probability for this hidden state
      p.hidden[h] = max(final.curr.prob)
    }
    
    # save the probabilities for this obswervation
    names(p.hidden)=hidden.states
    hidden.states.list[[i]]=p.hidden
  }
  
  # turn to a table and return
  hidden.states.viterbi.out = bind_rows(hidden.states.list)
  return(hidden.states.viterbi.out)
}

# run viterbi on our data
viterbi.tbl = My.Viterbi(
  emission.probabilities = emission.probabilities,
  transition.probabilities = transition.probabilities,
  priors = priors, 
  visible = visible
) 

# define a function to get optimal hidden state path
get.viterbi.path = function(
  viterbi.tbl=viterbi.tbl
)
{
  # get hidden state
  hidden.states = colnames(viterbi.tbl)
  
  # vector to save hidden state sequence
  viterbi.sequence=c()
  for (i in 1:dim(viterbi.tbl)[1]) 
  {
    # get probabilities for current row
    tmp=viterbi.tbl[i,]
    
    # get current state with maximum probability
    curr.state = hidden.states[tmp[1,]==max(tmp)]
    
    # resolve ties if two states are equally likely
    if(length(curr.state) > 1)
    {
      # pick randomly if ties
      curr.state=sample(
        x=curr.state,
        size = 1,
        replace = TRUE
      )
    }
    
    # add to sequence of hidden states
    viterbi.sequence=c(viterbi.sequence, curr.state)
  }
  # return the determined sequence
  return(viterbi.sequence)
}

# run
viterbi.path = get.viterbi.path(viterbi.tbl = viterbi.tbl)

# look at how many spots our code got correct
sum(data$Hidden == viterbi.path)

## [1] 19

This Viterbi algorithm is already implemented in the HMM R-Package. Let’s use the R-Package and compare to our implementation.

# set hmm object
hmm =initHMM(
    c("A", "B"),
    c(1, 2, 3), 
    startProbs = priors,
    transProbs = transition.probabilities,
    emissionProbs = t(emission.probabilities)
)

# package viterbi
true.viterbi = viterbi(
  hmm = hmm,
  data$Visible
)

# how many did the package get right
sum(true.viterbi == simulated.hmm$states)

## [1] 21

As we can see, our results are comparable to those delivered by packages. Differences can arise as two paths can be equally likely, and in our case, we are resolving ties randomly. We can compare further our algorithm implementation and the R-Package implementation by using simulation of HMM data, which can also be done with the HMM package. Let’s simulate a large sequence from the proposed HMM!

# simulate
simulated.hmm=simHMM(
  hmm = hmm,
  length = 10000
)

# We ca get inputs for our viterbi implementation from the HMM object
priors = hmm$startProbs
transition.probabilities = hmm$transProbs
emission.probabilities =  hmm$emissionProbs

# Our algorithm needs emission probabilities in correct format
emission.probabilities = data.frame(
  A = emission.probabilities['A',],
  B = emission.probabilities['B',]
)

# run my viterbi implementation
new.viterbi.tbl = My.Viterbi(
  emission.probabilities = emission.probabilities,
  transition.probabilities = transition.probabilities,
  priors = priors,
  visible =simulated.hmm$observation
)
new.viterbi.path = get.viterbi.path(viterbi.tbl = new.viterbi.tbl)

# package viterbi
true.viterbi = viterbi(
  hmm = hmm,
  simulated.hmm$observation
)

# compare to truth
mine=sum(new.viterbi.path == simulated.hmm$states)
pack=sum(true.viterbi == simulated.hmm$states)
mine

## [1] 5329

pack

## [1] 5237

# compare between implementations
sum(true.viterbi == new.viterbi.path)

## [1] 8332

As we can see, from data generated from the same Hidden Markov Model that is used as input to both implementations of the Viterbi algorithm, the R-Package and our implementation returns pretty similar results.

Finally, the inputs (priors, transitions, emissions) can be easily estimated from data

Example getting inputs from data

# get an example of training data from simulated data
training=data.frame(
  Hidden = simulated.hmm$states,
  Visible = simulated.hmm$observation
)
head(training)

##   Hidden Visible
## 1      B       1
## 2      A       3
## 3      A       2
## 4      A       3
## 5      B       2
## 6      B       3

# function to get probabilities from data
get.inputs.from.data = function(
  training = training
)
{
  # get priors 
  priors=table(training$Hidden)/dim(training)[1]
  priors
  
  # get transitions
  all.transitions=c()
  for (i in 1:(dim(training)[1]-1) ) 
  {
    # get both states
    first.state=training$Hidden[i]
    second.state=training$Hidden[i+1]
    
    # get pasted
    combination=paste(
      first.state,
      second.state,
      sep = '>'
    )
    all.transitions[i]=combination
  }
  transition.probabilities = table(all.transitions)/length(all.transitions)
  transition.probabilities
  
  # get transition probability matrix
  transition.probabilities  = matrix(
    data = transition.probabilities,
    nrow = 2,
    ncol = 2,
    byrow = TRUE
  )
  colnames(transition.probabilities)=c('A','B')
  rownames(transition.probabilities)=c('A','B')
  
  # emission probabilities
  emission.probabilities = table(training)/dim(training)[1]
  emission.probabilities
  
  # emission probabilities in correct format
  emission.probabilities = data.frame(
    A = emission.probabilities['A',],
    B = emission.probabilities['B',]
  )
  return.list=list(
    'emission.probabilities'=emission.probabilities,
    'transition.probabilities'=transition.probabilities,
    'priors'=priors
  )
  return(return.list)
}

# get inputs
inputs = get.inputs.from.data(training = training)

# generate some new data (test data)
simulated.hmm=simHMM(
  hmm = hmm,
  length = 10000
)

# run my viterbi
new.viterbi.tbl = My.Viterbi(
  emission.probabilities = inputs$emission.probabilities,
  transition.probabilities = inputs$transition.probabilities,
  priors = inputs$priors,
  visible = simulated.hmm$observation
)
new.viterbi.path = get.viterbi.path(viterbi.tbl = viterbi.tbl)

# package viterbi
true.viterbi = viterbi(
  hmm = hmm,
  observation = simulated.hmm$observation
)

# compare to truth
mine=sum(new.viterbi.path == simulated.hmm$states)
pack=sum(true.viterbi == simulated.hmm$states)
mine

## [1] 4969

pack

## [1] 5241

As expected, our implementation performed comparable to the R-Package, even though we estimated our probabilities from the data, and did not use the true HMM model as input (as we did for the R-Package). This is a good support that our estimations were good. Of course, the training data and the testing data were generated from the same HMM, which is an assumption that needs to also be true in real data analysis.