In jonotuke/HMMCPP: Implementation of HMM package in Rcpp

Installing lib

Use the command

devtools::install_github("jonotuke/HMMCPP")

Get libraries

library(HMMCPP)
library(tidyverse)

Crooked casino example

A crooked casino has two dice. One is fair and one is weighted so that 6 appears more often. The p.m.f for each is

| X | 1 | 2 | 3 | 4 | 5 | 6 | | --- | --- | --- | --- | --- | --- | --- | | Fair prob | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | 1/6 | | Biased prob | 1/10 | 1/10 | 1/10 | 1/10 | 1/10 | 1/2 |

When you enter the casino each die is equally likely to be used so

$$ P(F) = P(B) = 1/2, $$

where $F$ is the event that the fair die is used, and $B$ is the event that the biased die is used.

After each roll, there is a 10% chance that the die is changed giving the one-step transistion matrix:

| | F | B | | --- | --- | --- | | F | 0.9 | 0.1 | | B | 0.1 | 0.9 |

All we observe if the number rolled.

Set up in R

start_prob <- rep(1/2,2)
symbols <- 1:6
states <- c("F","B")
A <- matrix(c(0.9, 0.1, 0.1, 0.9), ncol = 2, byrow = TRUE)
A
B <- matrix(c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6,
              0.1, 0.1, 0.1, 0.1, 0.1, 0.5), 
            ncol = 6, byrow = TRUE)
B

Simulate some data

set.seed(2016)
crooked_casino <- sim_HMM(initial_prop = start_prob, 
                       n_sims = 500, 
                       A = A, 
                       state_sym = states, 
                       obs_sym = symbols, 
                       B = B)
crooked_casino

Look at the patterns

crooked_casino_df <- data_frame(state = crooked_casino$states, 
                             obs = crooked_casino$observation, 
                             time = seq_along(state))
crooked_casino_df

crooked_casino_df %>% 
  ggplot(aes(obs, fill = obs)) + 
  geom_bar(show.legend = FALSE) + 
  facet_wrap(~state)

Use package to get estimated states

estimated_states  <- Rcpp_viterbi(states = states, 
                                  symbols = symbols, 
                                  pi = start_prob, 
                                  A = A, 
                                  B = B, 
                                  obs = crooked_casino$observation)
table(estimated_states, crooked_casino$states)

Use package to get estimated parameters with Viterbi training

Rcpp_viterbi_training(states = states, 
                      symbols = symbols, 
                      pi = start_prob, 
                      A = A, 
                      B = B, 
                      obs = crooked_casino$observation)

Use package to get estimated parametes with Baum Welch

Rcpp_baum_welch(states = states, 
                symbols = symbols, 
                pi = start_prob, 
                A = A, 
                B = B, 
                obs = crooked_casino$observation)