A review of HiddenMarkov
The HiddenMarkov package needs some clarifica
I’m very pleased that Mr. David Harte has taken it upon himself to write an R package to perform computations for the Hidden Markov Model. Unfortunately, the documentation leaves me a bit confused about how to use his software under certain circumstances.
It is my goal to use the Viterbi algorithm to reproduce the state transitions of my Hidden Markov Model given the observations generated by the model. Specifically, I want to use my observation sequence to generate the state transition sequence given the observation sequence generated by a multinomial distribution in my Hidden Markov Model. This is known as the second of the three basic problems of a Hidden Markov Model.
To accomplish my goal of reproducing the state transition sequence given the generated observation sequence, I wrote this R script:
################################################################################
# observationCount is the number of observations we wish to generate.
observationCount <- 100
# observations is the observation vector. The valid values of observations are
# in V.
observations <- rep (LETTERS[24], observationCount)
numericObservations <- rep (0, observationCount)
# states is the actual state sequence used to generate the observation vector.
# The valid values of states are in {1, 2, ..., N}.
states <- rep (0, observationCount)
# V is a representation of the valid distinct observations in the model.
V <- c (LETTERS[c (18, 7, 2)]) # R, G, B
# delta is the set of initial probabilities for starting out in state S(1), S(2),
# ..., and S(N).
delta <- c (0.15, 0.25, 0.3, 0.3)
# N is the number of states in the model.
N <- length (delta)
# Pi is the state transition matrix. Each row represents the probabilities of
# transitioning from state S(i) to S(j), 1 <= i, j <= N.
Pi <- matrix (c (
0.2, 0.15, 0.1, 0.55,
0.25, 0.25, 0.25, 0.25,
0.2, 0.3, 0.3, 0.2,
0.25, 0.15, 0.5, 0.1),
nrow = N, byrow = T)
# The distinct observations are listed in V. V is of size M, and the number of
# columns in B is M. B is the matrix of multinomially distributed observation
# probabilities in state S(i), 1 <= i <= N.
B <- matrix (c (
10/100, 40/100, 50/100,
45/90, 10/90, 35/90,
12/112, 61/112, 39/112,
41/108, 23/108, 44/108),
nrow = N, byrow = T)
# M is the number of possible distinct observations in the experiment. V is of
# size M.
M <- ncol (B)
################################################################################
firstState <- function (P) {
# P: The vector of initial state probabilities.
# We select the initial state.
cumulativeP <- c (P[1], sum (P[1:2]), sum (P[1:3]), sum (P[1:4]))
x <- runif (1)
if (x < cumulativeP[1])
{
s <- 1
}
else if (x < cumulativeP[2])
{
s <- 2
}
else if (x < cumulativeP[3])
{
s <- 3
}
else
{
s <- 4
}
return (s)
}
################################################################################
nextState <- function (s, P) {
# s: The current state. This state maps into a row of the state transition
# matrix.
# P: The state transition matrix, i.e., the probability of transitioning from
# state i to state j, 1 <= i, j <= N.
# We select the next state given that we are currently in state s.
cumulativeP <- c (P[s, 1], sum (P[s, 1:2]), sum (P[s, 1:3]), sum (P[s, 1:4]))
q <- 0
x <- runif (1)
if (x <= cumulativeP[1])
{
q <- 1
}
else if (x <= cumulativeP[2])
{
q <- 2
}
else if (x <= cumulativeP[3])
{
q <- 3
}
else
{
q <- 4
}
return (q)
}
################################################################################
observation <- function (r, Obs) {
# r: The row of Obs which represents the current state, 1 <= r <= N.
# Obs: The N by M matrix of observation probabilities.
cumulativeObs <- c (Obs[r, 1], sum (Obs[r, 1:2]), sum (Obs[r, 1:3]))
obs <- LETTERS[24]
x <- runif (1)
if (x < cumulativeObs[1])
{
obs <- V[1]
}
else if (x < cumulativeObs[2])
{
obs <- V[2]
}
else
{
obs <- V[3]
}
return (obs)
}
################################################################################
convertToNumericObservations <- function (obs) {
# obs: The observation vector which is of size observationCount.
nObs <- rep (0, length (obs))
for (i in 1:length (obs)) {
if (obs[i] == V[1]) {
nObs[i] <- 1
}
else if (obs[i] == V[2]) {
nObs[i] <- 2
}
else {
nObs[i] <- 3
}
}
return (nObs)
}
################################################################################
states[1] <- firstState (delta)
for (i in 1:observationCount) {
observations[i] <- observation (states[i], B)
if (i < observationCount) {
states[i + 1] <- nextState (states[i], Pi)
}
}
numericObservations <- convertToNumericObservations (observations)
print (states)
print (observations)
print (numericObservations)
################################################################################
#
library (HiddenMarkov)
HMM <- dthmm (x = numericObservations, Pi = Pi, delta = delta, distn = "multinom",
pm = list (size = NULL), pn = list (prob = B), discrete = TRUE)
computedStates <- Viterbi (HMM)
################################################################################
All of this works just fine until I get to the final line: computedStates <- Viterbi (HMM) where I get the error message: "Error in if (any(prob < 0) || (s <- sum(prob)) == 0) stop("probabilities cannot be negative nor all 0.") : missing value where TRUE/FALSE needed"
I used the R debugger to trace through the code. Evidently, the observation matrix, B, is cycled through once. At the end of the matrix it grabs the next probability value which is outside the bounds of B. An invalid value is returned.
So my question is this: How can I use the HiddenMarkov package to reproduce the state transition vector? I’m sure that there is something that I’m unaware of that will allow me to use the Viterbi algorithm in this package.
Other than the minor weakness in the documentation, the HiddenMarkov package works very well.
