# Hidden Markov Models

## Simulating HMM

``````hidden_states = c("hot", "cold")

N = length(hidden_states)

# transition matrix of hidden states
A = rbind(
c(0.7, 0.3),
c(0.4, 0.6)
)

# initial distribution of hidden state
pi_0 = c(0.571, 0.429)

observed_states = c("small", "medium", "large")

M = length(observed_states)

# emission matrix: hidden -> observed
B = rbind(
c(0.1, 0.4, 0.5),
c(0.7, 0.2, 0.1)
)``````
``````T = 100

# Note: +1's below because first entries represent initial time 0
X = rep(NA, T + 1)
O = rep(NA, T + 1)

# initial state
X[1] = sample(1:N, 1, prob = pi_0)
O[1] = sample(1:M, 1, prob = B[X[1], ])

for (t in 1:T) {
# simulate current hidden state, given previous hidden state
X[t + 1] = sample(1:N, 1, prob = A[X[t], ])

# simulate current observed state, given current hidden state
O[t + 1] = sample(1:M, 1, prob = B[X[t + 1], ])

}

data.frame(t = 0:T, X_t = hidden_states[X], O_t = observed_states[O]) |>
t X_t O_t
0 hot large
1 hot large
2 hot large
3 hot large
4 hot large
5 cold large
6 cold medium
7 hot medium
8 hot medium
9 hot small
``````plot(0:T, X, type = "o", xlab = "t", ylab = "X_t", col = "blue",
main = "Hidden states (1 = hot, 2 = cold)")
plot(0:T, O, type = "o", xlab = "t", ylab = "O_t", col = "orange",
main = "Observed states (1 = small, 2 = medium, 3 = large)")``````

## Stationary/Marginal at single time

``````# hidden states
pi_sd = compute_stationary_distribution(A)

pi_sd``````
``````          [,1]      [,2]
[1,] 0.5714286 0.4285714``````
``````# observed states
pi_sd %*% B``````
``````          [,1]      [,2]      [,3]
[1,] 0.3571429 0.3142857 0.3285714``````

## HMM package

``````library(HMM)

# Initialise HMM
hmm = initHMM(States = hidden_states,
Symbols = observed_states,
transProbs = A,
emissionProbs = B,
startProbs = pi_sd)

hmm``````
``````\$States
[1] "hot"  "cold"

\$Symbols
[1] "small"  "medium" "large"

\$startProbs
hot      cold
0.5714286 0.4285714

\$transProbs
to
from   hot cold
hot  0.7  0.3
cold 0.4  0.6

\$emissionProbs
symbols
states small medium large
hot    0.1    0.4   0.5
cold   0.7    0.2   0.1``````

## HMM package simulation

``simHMM(hmm, 10)``
``````\$states
[1] "hot"  "cold" "hot"  "cold" "hot"  "cold" "cold" "cold" "cold" "hot"

\$observation
[1] "medium" "small"  "large"  "small"  "medium" "small"  "small"  "small"
[9] "small"  "medium"``````

## Likelihood of observations

``````# Sequence of observations
observations = c("small", "medium")

log_forward_probability = forward(hmm, observations)

forward_probability = exp(log_forward_probability)

forward_probability``````
``````      index
states          1          2
hot  0.05714286 0.06400000
cold 0.30000000 0.03942857``````
``colSums(forward_probability)``
``````        1         2
0.3571429 0.1034286 ``````
``````# Sequence of observations
observations = simHMM(hmm, 10)\$observation
observations``````
`````` [1] "large"  "large"  "medium" "small"  "medium" "small"  "large"  "medium"
[9] "small"  "medium"``````
``````log_forward_probability = forward(hmm, observations)

forward_probability = exp(log_forward_probability)

forward_probability``````
``````      index
states          1          2           3           4           5            6
hot  0.28571429 0.10857143 0.032182857 0.002566857 0.002327680 0.0002173669
cold 0.04285714 0.01114286 0.007851429 0.010056000 0.001360731 0.0010603200
index
states            7            8            9           10
hot  2.881424e-04 9.190230e-05 7.461376e-06 6.904504e-06
cold 7.014021e-05 2.570537e-05 3.009574e-05 4.059171e-06``````
``colSums(forward_probability)``
``````           1            2            3            4            5            6
3.285714e-01 1.197143e-01 4.003429e-02 1.262286e-02 3.688411e-03 1.277687e-03
7            8            9           10
3.582826e-04 1.176077e-04 3.755711e-05 1.096367e-05 ``````

## Decoding

``````observations = c("small", "medium")

viterbi(hmm, observations)``````
``[1] "cold" "hot" ``
``````# Sequence of observation
observations = simHMM(hmm, 100)\$observation

observations``````
``````  [1] "large"  "medium" "large"  "medium" "small"  "large"  "large"  "small"
[9] "small"  "small"  "medium" "medium" "large"  "large"  "small"  "medium"
[17] "medium" "large"  "medium" "small"  "large"  "medium" "medium" "small"
[25] "large"  "small"  "medium" "medium" "large"  "large"  "large"  "small"
[33] "large"  "small"  "small"  "small"  "medium" "large"  "large"  "large"
[41] "large"  "large"  "medium" "large"  "medium" "medium" "small"  "small"
[49] "large"  "small"  "small"  "medium" "large"  "small"  "medium" "small"
[57] "large"  "medium" "medium" "large"  "small"  "medium" "medium" "large"
[65] "large"  "medium" "large"  "small"  "small"  "large"  "small"  "medium"
[73] "large"  "small"  "medium" "small"  "medium" "medium" "medium" "medium"
[81] "small"  "small"  "small"  "medium" "medium" "small"  "small"  "large"
[89] "large"  "large"  "medium" "small"  "large"  "small"  "small"  "medium"
[97] "small"  "large"  "medium" "large" ``````
``viterbi(hmm, observations)``
``````  [1] "hot"  "hot"  "hot"  "hot"  "cold" "hot"  "hot"  "cold" "cold" "cold"
[11] "hot"  "hot"  "hot"  "hot"  "cold" "hot"  "hot"  "hot"  "hot"  "cold"
[21] "hot"  "hot"  "hot"  "cold" "hot"  "cold" "hot"  "hot"  "hot"  "hot"
[31] "hot"  "cold" "hot"  "cold" "cold" "cold" "hot"  "hot"  "hot"  "hot"
[41] "hot"  "hot"  "hot"  "hot"  "hot"  "hot"  "cold" "cold" "hot"  "cold"
[51] "cold" "hot"  "hot"  "cold" "cold" "cold" "hot"  "hot"  "hot"  "hot"
[61] "cold" "hot"  "hot"  "hot"  "hot"  "hot"  "hot"  "cold" "cold" "hot"
[71] "cold" "hot"  "hot"  "cold" "cold" "cold" "hot"  "hot"  "hot"  "hot"
[81] "cold" "cold" "cold" "hot"  "hot"  "cold" "cold" "hot"  "hot"  "hot"
[91] "hot"  "cold" "hot"  "cold" "cold" "cold" "cold" "hot"  "hot"  "hot" ``````

## Estimation

``````# Sequence of observation
observation = simHMM(hmm, 10000)\$observation

table(observation) / 10000``````
``````observation
large medium  small
0.3205 0.3211 0.3584 ``````
``````# Baum-Welch
baumWelch(hmm, observation, maxIterations = 10)``````
``````\$hmm
\$hmm\$States
[1] "hot"  "cold"

\$hmm\$Symbols
[1] "small"  "medium" "large"

\$hmm\$startProbs
hot      cold
0.5714286 0.4285714

\$hmm\$transProbs
to
from         hot      cold
hot  0.6990122 0.3009878
cold 0.4018597 0.5981403

\$hmm\$emissionProbs
symbols
states      small    medium     large
hot  0.09936791 0.4157953 0.4848368
cold 0.70413347 0.1947090 0.1011575

\$difference
[1] 0.018781129 0.002290127 0.002154031 0.002030550 0.001917701 0.001813952
[7] 0.001718117 0.001629262 0.001546632 0.001469609``````

## All or nothing

``````library(HMM)

hidden_states = c("learning", "mastery")

observed_states = c("incorrect", "correct")

A = rbind(
c(0.8, 0.2),
c(0.001, 0.999)
)

B = rbind(
c(0.5, 0.5),
c(0.1, 0.9)
)

pi_0 = c(0.999, 0.001)

# Initialise HMM
hmm = initHMM(States = hidden_states,
Symbols = observed_states,
transProbs = A,
emissionProbs = B,
startProbs = pi_0)

hmm``````
``````\$States
[1] "learning" "mastery"

\$Symbols
[1] "incorrect" "correct"

\$startProbs
learning  mastery
0.999    0.001

\$transProbs
to
from       learning mastery
learning    0.800   0.200
mastery     0.001   0.999

\$emissionProbs
symbols
states     incorrect correct
learning       0.5     0.5
mastery        0.1     0.9``````
``simHMM(hmm, 10)``
``````\$states
[1] "learning" "learning" "mastery"  "mastery"  "mastery"  "mastery"
[7] "mastery"  "mastery"  "mastery"  "mastery"

\$observation
[1] "correct" "correct" "correct" "correct" "correct" "correct" "correct"
[8] "correct" "correct" "correct"``````
``````# Sequence of observation
observations = simHMM(hmm, 100)\$observation

observations``````
``````  [1] "incorrect" "correct"   "correct"   "correct"   "correct"   "correct"
[7] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[13] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[19] "correct"   "incorrect" "correct"   "correct"   "correct"   "correct"
[25] "correct"   "correct"   "correct"   "incorrect" "correct"   "correct"
[31] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[37] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[43] "correct"   "correct"   "correct"   "incorrect" "correct"   "correct"
[49] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[55] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[61] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[67] "correct"   "correct"   "correct"   "correct"   "correct"   "incorrect"
[73] "correct"   "incorrect" "correct"   "correct"   "incorrect" "correct"
[79] "incorrect" "correct"   "correct"   "correct"   "correct"   "incorrect"
[85] "incorrect" "correct"   "correct"   "incorrect" "incorrect" "incorrect"
[91] "correct"   "correct"   "correct"   "correct"   "correct"   "correct"
[97] "correct"   "correct"   "correct"   "correct"  ``````
``viterbi(hmm, observations)``
``````  [1] "learning" "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[7] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[13] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[19] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[25] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[31] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[37] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[43] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[49] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[55] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[61] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[67] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[73] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[79] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[85] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[91] "mastery"  "mastery"  "mastery"  "mastery"  "mastery"  "mastery"
[97] "mastery"  "mastery"  "mastery"  "mastery" ``````
``````# Sequence of observation
observation = simHMM(hmm, 10000)\$observation
table(observation) / 10000``````
``````observation
correct incorrect
0.8988    0.1012 ``````
``````# Baum-Welch
baumWelch(hmm, observation, maxIterations = 10)``````
``````\$hmm
\$hmm\$States
[1] "learning" "mastery"

\$hmm\$Symbols
[1] "incorrect" "correct"

\$hmm\$startProbs
learning  mastery
0.999    0.001

\$hmm\$transProbs
to
from          learning   mastery
learning 0.895489932 0.1045101
mastery  0.001329809 0.9986702

\$hmm\$emissionProbs
symbols
states      incorrect   correct
learning 0.45947571 0.5405243
mastery  0.09629494 0.9037051

\$difference
[1] 0.112359065 0.040630020 0.032319113 0.024353848 0.018073108 0.013406443
[7] 0.010007670 0.007541187 0.005744390 0.004425991``````