A Revealing Introduction to Hidden Markov Models

Mark Stamp
January 18, 2004
1 A simple example
Suppose we want to determine the average annual temperature at a particular location on
earth over a series of years. To make it interesting, suppose the years we are concerned with
lie in the distant past, before thermometers were invented. Since we cant go back in time,
we instead look for indirect evidence of the temperature.
To simplify the problem, we only consider two annual temperatures, hot and cold.
Suppose that modern evidence indicates that the probability of a hot year followed by another
hot year is 0.7 and the probability that a cold year is followed by another cold year is 0.6.
Well assume that these probabilities held in the distant past as well. The information so
far can be summarized as
H C
H
C
_
0.7 0.3
0.4 0.6
_
(1)
where H is hot and C is cold.
Also suppose that current research indicates a correlation between the size of tree growth
rings and temperature. For simplicity, we only consider three dierent tree ring sizes, small,
medium and large, or S, M and L. Conceivably, the probabilistic relationship between
temperature and tree ring sizes could be given by
S M L
H
C
_
0.1 0.4 0.5
0.7 0.2 0.1
_
.
(2)
For this system, the state is the average annual temperatureeither H or C. The tran-
sition from one state to the next is a Markov process (of order one), since the next state
depends only on the current state and the xed probabilities in (1). However, the actual
states are hidden since we cant directly observe the temperature in the past.
Although we cant observe the state (temperature) we can observe the size of tree rings.
From (2), tree rings provide us with probabilistic information regarding the temperature.
Since the states are hidden, this type of system is known as a Hidden Markov Model (HMM).
1
Our goal is to make eective and ecient use of the observable information in order to gain
insight into various aspects of the Markov process.
The state transition matrix
A =
_
0.7 0.3
0.4 0.6
_
(3)
comes from (1) and the observation matrix
B =
_
0.1 0.4 0.5
0.7 0.2 0.1
_
. (4)
is from (2). Perhaps there is additional evidence that the initial state distribution, denoted
by , is
=
_
0.6 0.4
_
. (5)
The matrices, , A and B are row stochastic, meaning that each element is a probability
and the elements of each row sum to 1.
Now consider a particular four-year period of interest where we observe the series of tree
rings S, M, S, L. Letting 0 represent S, 1 represent M and 2 represent L, this observation
sequence is
O = (0, 1, 0, 2). (6)
We might want to determine the most likely state sequence of the Markov process given
the observations (6). This is not quite as clear-cut as it seems, since there are dierent
possible interpretations of most likely. On the one hand, we could reasonably dene
most likely as the state sequence with the highest probability from among all possible state
sequences of length four. Dynamic programming (DP) eciently nds this particular most
likely solution. On the other hand, we might reasonably dene most likely as the state
sequence that maximizes the expected number of correct states. HMMs nd this most likely
sequence.
The DP solution and the HMM solution are not necessarily the same. For example, the
DP solution must have valid state transitions, while this is not necessarily the case for the
HMMs. And even if all state transitions are valid, the HMM solution can still dier from
the DP solution (as we illustrate in an example below).
Next, we present one of the most challenging aspects of HMMsthe notation. Then
we discuss the three fundamental problems related to HMMs and give algorithms for their
ecient solutions. We also consider some critical computational issue that must be addressed
when writing any HMM computer program. We conclude with a substantial example that
does not require any specialized knowledge, yet nicely illustrates the strength of the HMM
approach. Rabiner [1] is the best source for further introductory information on HMMs.
2
2 Notation
Let
T = the length of the observation sequence
N = the number of states in the model
M = the number of observation symbols
Q = {q
0
, q
1
, . . . , q
N1
} = the states of the Markov process
V = {0, 1, . . . , M 1} = set of possible observations
A = the state transition probabilities
B = the observation probability matrix
= the initial state distribution
O = (O
0
, O
1
, . . . , O
T1
) = observation sequence.
The observations are always denoted by {0, 1, . . . , M 1}, since this simplies the notation
with no loss in generality. Then O
i
V for i = 0, 1, . . . , T 1.
A generic hidden Markov model is illustrated in Figure 1, where the X
i
are the hidden
states and all other notation is as given above. The Markov processwhich is hidden behind
the dashed lineis determined by the initial state X
0
and the A matrix. We are only able
to observe the O
i
, which are related to the actual states of the Markov process by the
matrices B and A.
Markov process:
X
0
X
1
X
2
X
T1
-
A
-
A
-
A
-
A
?
B
?
B
?
B
?
B
Observations:
O
0
O
1
O
2
O
T1
Figure 1: Hidden Markov Model
For the temperature example of the previous sectionwith the observations sequence
given in (6)we have T = 4, N = 2, M = 3, Q = {H, C}, V = {0, 1, 2} (where we let 0, 1, 2
represent small, medium and large tree rings). In this case, the matrices A, B and
are given by (3), (4) and (5), respectively.
The matrix A = {a
ij
} is N N with
a
ij
= P(state q
j
at t + 1 | state q
i
at t)
and A is row stochastic. Note that the probabilities a
ij
are independent of t. The matrix
B = {b
j
(k)} is an N M with
b
j
(k) = P(observation k at t | state q
j
at t).
3
As with A, the matrix B is row stochastic and the probabilities b
j
(k) are independent of t.
The unusual notation b
j
(k) is standard in the HMM world.
An HMM is dened by A, B and (and, implicitly, by the dimensions N and M). The
HMM is denoted by = (A, B, ).
Consider a state generic sequence of length four
X = (x
0
, x
1
, x
2
, x
3
)
with corresponding observations
O = (O
0
, O
1
, O
2
, O
3
).
Then
x
0
is the probability of starting in state x
0
. Also, b
x
0
(O
0
) is the probability of initially
observing O
0
and a
x
0
,x
1
is the probability of transiting from state x
0
to state x
1
. Continuing,
we see that the probability of the state sequence X is given by
P(X) =
x
0
b
x
0
(O
0
)a
x
0
,x
1
b
x
1
(O
1
)a
x
1
,x
2
b
x
2
(O
2
)a
x
2
,x
3
b
x
3
(O
3
). (7)
Consider the temperature example in Section 1 with observation sequence O = (0, 1, 0, 2),
as given in (6). Using (7) we can compute, say,
P(HHCC) = 0.6(0.1)(0.7)(0.4)(0.3)(0.7)(0.6)(0.1) = 0.000212.
Similarly, we can compute the probability of each of the possible state sequences of length
four, assuming the xed observation sequence (6). We have listed these results in Table 1,
where the probabilities in the last column are normalized so that they sum to 1.
To nd the optimal state sequence in the dynamic programming (DP) sense, we simply
choose the sequence with the highest probability, namely, CCCH. To nd the optimal
sequence in the HMM sense, we choose the most probable symbol at each position. To this
end we sum the probabilities in Table 1 that have an H in the rst position. Doing so, we nd
the (normalized) probability of H in the rst position is 0.18817 and hence the probability
of C in the rst position is 0.81183. The HMM therefore chooses the rst element of the
optimal sequence to be C. We repeat this for each element of the sequence, obtaining the
probabilities in Table 2.
From Table 2 we nd that the optimal sequencein the HMM senseis CHCH. Note
that in this example, the optimal DP sequence diers from the optimal HMM sequence, even
though all state transitions are valid.
3 The three problem
3.1 Problem 1
Given the model = (A, B, ) and a sequence of observations O, nd P(O| ). Here, we
want to determine the likelihood of the observed sequence O, given the model.
4
normalized
state probability probability
HHHH .000412 .042743
HHHC .000035 .003664
HHCH .000706 .073274
HHCC .000212 .021982
HCHH .000050 .005234
HCHC .000004 .000449
HCCH .000302 .031403
HCCC .000091 .009421
CHHH .001098 .113982
CHHC .000094 .009770
CHCH .001882 .195398
CHCC .000564 .058619
CCHH .000470 .048849
CCHC .000040 .004187
CCCH .002822 .293096
CCCC .000847 .087929
Table 1: State sequence probabilities
element
0 1 2 3
P(H) 0.188170 0.519432 0.228878 0.803979
P(C) 0.811830 0.480568 0.771122 0.196021
Table 2: HMM probabilities
3.2 Problem 2
Given = (A, B, ) and an observation sequence O, nd an optimal state sequence for the
underlying Markov process. In other words, we want to uncover the hidden part of the
Hidden Markov Model. This type of problem is discussed in Section 1.
3.3 Problem 3
Given an observation sequence O and the dimensions N and M, nd the model = (A, B, )
that maximizes the probability of O. This can be viewed as training the model to best t
the observed data. Alternatively, we can view this as a (discrete) hill climb on the parameter
space represented by A, B and .
5
3.4 Discussion
Consider speech recognitionwhich happens to be one of the best-known applications of
HMMs. We can use the solution to Problem 3 to train an HMM, say,
0
to recognize the
spoken word no and train another HMM, say,
1
to recognize the spoken word yes.
Then given an unknown spoken word, we can use the solution to Problem 1 to score this
word against
0
and also
1
to determine whether it is more likely no, yes or neither. In
this case, we dont need to solve Problem 2, but it is possible that such a solutionwhich
uncovers the hidden statesmight provide additional insight into the underlying speech
model.
4 The three solutions
4.1 Solution to Problem 1
Let = (A, B, ) be a given model and let O = (O
0
, O
1
, . . . , O
T1
) be a series of observations.
We want to nd P(O| ).
Let X = (x
0
, x
1
, . . . , x
T1
) be a state sequence. Then by the denition of B we have
P(O| X, ) = b
x
0
(O
0
)b
x
1
(O
1
) b
x
T1
(O
T1
)
and by the denition of and A it follows that
P(X | ) =
x
0
a
x
0
,x
1
a
x
1
,x
2
a
x
T2
,x
T1
.
Since
P(O, X | ) =
P(O X )
P()
and
P(O| X, )P(X | ) =
P(O X )
P(X )

P(X )
P()
=
P(O X )
P()
we have
P(O, X | ) = P(O| X, )P(X | ).
By summing over all possible state sequences we obtain
P(O| ) =

X
P(O, X | )
=

X
P(O| X, )P(X | )
=

x
0
b
x
0
(O
0
)a
x
0
,x
1
b
x
1
(O
1
) a
x
T2
,x
T1
b
x
T1
(O
T1
).
However, this direct computation is generally infeasible, since it requires about 2TN
T
multi-
plications. The strength of the HMM approach derives largely from the fact that there exist
an ecient algorithm to achieve the same result.
6
To nd P(O| ), the so-called forward algorithm, or -pass, is used. For t = 0, 1, . . . , T1
and i = 0, 1, . . . , N 1, dene

t
(i) = P(O
0
, O
1
, . . . , O
t
, x
t
= q
i
| ). (8)
Then
t
(i) is the probability of the partial observation sequence up to time t, where the
underlying Markov process in state q
i
at time t.
The key result is that the
t
(i) can be computed recursively as
1. Let
0
(i) =
i
b
i
(O
0
), for i = 0, 1, . . . , N 1
2. For t = 1, 2, . . . , T 1 and i = 0, 1, . . . , N 1, compute

t
(i) =
_
_
N1

j=0

t1
(j)a
ji
_
_
b
i
(O
t
)
3. Then from (8) it is clear that
P(O| ) =
N1

i=0

T1
(i).
The -pass computation requires N
2
T multiplications, as opposed to more than 2TN
T
for
the nave approach.
4.2 Solution to Problem 2
Given the model = (A, B, ) and a sequence of observations O, our goal is to nd the
most likely state sequence. As mentioned above, there are dierent possible interpretations
of most likely. For HMMs we maximize the expected number of correct states. In contrast,
a dynamic program nds the highest scoring overall path. As we have seen, these solutions
are not necessarily the same.
First, we dene the backward algorithm, or -pass. This is analogous to the -pass
discussed above, except that it starts at the end and works back toward the beginning.
For t = 0, 1, . . . , T 1 and i = 0, 1, . . . , N 1, dene

t
(i) = P(O
t+1
, O
t+2
, . . . , O
T1
| x
t
= q
i
, ).
Then the
t
(i) can be computed recursively as
1. Let
T1
(i) = 1, for i = 0, 1, . . . , N 1.
2. For t = T 2, T 1, . . . , 0 and i = 0, 1, . . . , N 1 compute

t
(i) =
N1

j=0
a
ij
b
j
(O
t+1
)
t+1
(j).
7
For t = 0, 1, . . . , T 2 and i = 0, 1, . . . , N 1, dene

t
(i) = P(x
t
= q
i
| O, ).
Since
t
(i) measures the relevant probability up to time t and
t
(i) measures the relevant
probability after time t,

t
(i) =

t
(i)
t
(i)
P(O| )
.
Recall that the denominator P(O| ) is obtained by summing
T1
(i) over i. From the
denition of
t
(i) it follows that the most likely state at time t is the state q
i
for which
t
(i)
is maximumwhere the maximum is taken over the index i.
4.3 Solution to Problem 3
Here we want to adjust the model parameters to best t the observations. The sizes of the
matrices (N and M) are xed but the elements of A, B and are free, subject only to the
row stochastic condition. The fact that we can re-estimate the model itself is one of the
more amazing aspects of HMMs.
For t = 0, 1, . . . , T 2 and i, j {0, 1, . . . , N 1}, dene di-gammas as

t
(i, j) = P(x
t
= q
i
, x
t+1
= q
j
| O, ).
Then
t
(i, j) is the probability of being in state q
i
at time t and transiting to state q
j
at
time t + 1. The di-gammas can be written in terms of , , A and B as

t
(i, j) =

t
(i)a
ij
b
j
(O
t+1
)
t+1
(j)
P(O| )
.
The
t
(i) and
t
(i, j) (or di-gamma) are related by

t
(i) =
N1

j=0

t
(i, j).
Given the and di-gamma we verify below that the model = (A, B, ) can be re-
estimated as
1. For i = 0, 1, . . . , N 1, let

i
=
0
(i) (9)
2. For i = 0, 1, . . . , N 1 and j = 0, 1, . . . , N 1, compute
a
ij
=
T2

t=0

t
(i, j)
_
T2

t=0

t
(i) . (10)
8
3. For j = 0, 1, . . . , N 1 and k = 0, 1, . . . , M 1, compute
b
j
(k) =

t{0,1,...,T2}
O
t
=k

t
(j)
_
T2

t=0

t
(j) . (11)
The numerator of the re-estimated a
ij
can be seen to give the expected number of tran-
sitions from state q
i
to state q
j
, while the denominator is the expected number of transitions
from q
i
to any state. Then the ratio is the probability of transiting from state q
i
to state q
j
,
which is the desired value of a
ij
.
The numerator of the re-estimated b
j
(k) is the expected number of times the model is
in state q
j
with observation k, while the denominator is the expected number of times the
model is in state q
j
. The ratio is the probability of observing symbol k, given that the model
is in state q
j
, which is the desired value of b
j
(k).
Re-estimation is an iterative process. First, we initialize = (A, B, ) with a best
guess or, if no reasonable guess is available, we choose random values such that
i
1/N
and a
ij
1/N and b
j
(k) 1/M. Its critical that A, B and be randomized, since exactly
uniform values result in a local maximum from which the model cannot climb. As always, ,
A and B must be row stochastic.
The process can be summarized as
1. Initialize, = (A, B, ).
2. Compute
t
(i),
t
(i),
t
(i, j) and
t
(i).
3. Re-estimate the model = (A, B, ).
4. If P(O| ) increases, goto 2.
Of course, it might be desirable to stop if P(O| ) does not increase by at least some
predetermined threshold and/or to set a maximum number of iterations.
5 Dynamic programming
Before completing our discussion of the elementary aspects of HMMs, we make a brief detour
to show the relationship between dynamic programming (DP) and HMMs. In fact, a DP
can be viewed as an -pass where sum is replaced by max. More precisely, for , A
and B as above, the dynamic programming algorithm can be stated as
1. Let
0
(i) =
i
b
i
(O
0
), for i = 0, 1, . . . , N 1.
2. For t = 1, 2, . . . , T 1 and i = 0, 1, . . . , N 1, compute

t
(i) = max
j{0,1,...,N1}
[
t1
(j)a
ji
b
i
(O
t
)]
9
At each successive t, the DP determines the probability of the best path ending at each of
the states i = 0, 1, . . . , N 1. Consequently, the probability of the best overall path is
max
j{0,1,...,N1}
[
T1
(j)]. (12)
Be sure to note that (12) only gives the optimal probability, not the optimal path itself.
By keeping track of each preceding state, the DP procedure can be used to recover the
optimal path by tracing back from the highest-scoring nal state.
For example, consider the example of Section 1. The initial probabilites are
P(H) =
0
b
0
(0) = 0.6(0.1) = 0.06 and P(T) =
1
b
1
(0) = 0.4(0.7) = 0.28.
The probabilities of each path of length two are
P(HH) = 0.06(0.7)(0.4) = 0.0168
P(HC) = 0.06(0.3)(0.2) = 0.0036
P(CH) = 0.28(0.4)(0.4) = 0.0448
P(CC) = 0.28(0.6)(0.2) = 0.0336
and hence the best (most probable) path of length two ending with H is CH while the
best path of length two ending with C is CC. Continuing, we construct the the diagram
in Figure 2 one level at a time, where each arrow points to the preceding element in the
optimal path up to that particular state. Note that at each stage, the dynamic programming
algorithm only needs to maintain the highest-scoring paths at each possible statenot a list
of all possible paths. This is the key to the eciency of the DP algorithm.
H
.06
C
.28
H
.0448

=
C
.0336

H
.003136

C
.014112

H
.002822

=
C
.000847

Figure 2: Dynamic programming

0
(i) = log[
i
(O
0
)], for i = 0, 1, . . . , N 1.
2. For t = 1, 2, . . . , T 1 and i = 0, 1, . . . , N 1, compute

t
(i) = max
j{0,1,...,N1}
_

t1
(j) + log[a
ji
] + log[b
i
(O
t
)]
_
.
10
In this case, the optimal score is
max
j{0,1,...,N1}
[

T1
(j)].
Of course, additional bookkeeping is still required in order nd the optimal path.
6 HMM Scaling
The three HMM solutions in Section 4 all require computations involving products of prob-
abilities. It is easy to see, for example, that
t
(i) tends to 0 exponentially as T increases.
Therefore, any attempt to implement the formulae as given above will inevitably result in
underow. The solution to this underow problem is to scale the numbers. However, care
must be taken to insure that, for example, the re-estimation formulae remain valid.
First, consider the computation of
t
(i). The basic recurrence is

t
(i) =
N1

j=0

t1
(j)a
ji
b
i
(O
t
).
It seems sensible to normalize each
t
(i) by dividing by the sum (over j) of
t
(j). However,
we must verify that the re-estimation formulae hold.
For t = 0, let
0
(i) =
0
(i) for i = 0, 1, . . . , N 1. Then let c
0
= 1/

N1
j=0

0
(j) and,
nally,
0
(i) = c
0

0
(i) for i = 0, 1, . . . , N 1. Then for each t = 1, 2, . . . , T 1 do the
following.
1. For i = 0, 1, . . . , N 1, compute

t
(i) =
N1

j=0

t1
(j)a
ji
b
i
(O
t
).
2. Let
c
t
=
1
N1

j=0

t
(j)
.
3. For i = 0, 1, . . . , N 1, compute

t
(i) = c
t

t
(i).
Clearly,
0
(i) = c
0

0
(i). Suppose that

t
(i) = c
0
c
1
c
t

t
(i). (13)
11
Then

t+1
(i) = c
t+1

t+1
(i)
= c
t+1
N1

j=0

t
(j)a
ji
b
i
(O
t+1
)
= c
0
c
1
c
t
c
t+1
N1

j=0

t
(j)a
ji
b
i
(O
t+1
)
= c
0
c
1
c
t+1

t+1
(i)
and hence (13) holds, by induction, for all t.
From (13) and the denitions of and it follows that

t
(i) =

t
(i)
N1

j=0

t
(j)
(14)
and hence
t
(i) are the desired scaled values of
t
(i) for all t.
As a consequence of (14),
N1

j=0

T1
(j) = 1.
Also, from (13) we have
N1

j=0

T1
(j) = c
0
c
1
c
T1
N1

j=0

T1
(j)
= c
0
c
1
c
T1
P(O| ).
Combining these results gives
P(O| ) =
1
T1

j=0
c
j
.
To avoid underow, we instead compute
log[P(O| )] =
T1

j=0
log c
j
. (15)
The same scale factor is used for
t
(i) as was used for
t
(i), namely c
t
, so that

t
(i) =
c
t

t
(i). We then compute
t
(i, j) and
t
(i) using the formulae of the previous section
with
t
(i) and

t
(i) in place of
t
(i) and
t
(i), respectively. These values are then used
to re-estimate , A and B.
By writing the original re-estimation formulae (9) and (10) and (11) directly in terms
of
t
(i) and
t
(i), it is an easy exercise to show that the re-estimated and A and B
12
are exact when
t
(i) and

t
(i) are used in place of
t
(i) and
t
(i). Furthermore, P(O| )
isnt required in the re-estimation formulae, since in each case it cancels in the numerator
and denominator. Therefore, (15) can be used to verify that P(O| ) is increasing at each
iteration. Fortunately, we have no need for the actual value of P(O| ), the calculation of
which would inevitably result in underow.
7 Putting it all together
Here we give complete pseudo-code for solving Problem 3, including scaling. This pseudo-
code also provides everything needed to solve Problems 1 and 2.
The values N and M are xed and the T observations
O = (O
0
, O
1
, . . . , O
T1
)
are assumed known.
1. Initialization:
Select initial values for the matrices A, B and , where is 1 N, while A = {a
ij
}
is N N and B = {b
j
(k)} is N M, and all three matrices are row-stochastic. If
known, use reasonable approximations for the matrix values, otherwise let
i
1/N
and a
ij
1/N and b
j
(k) 1/M. Be sure that each row sums to 1 and the elements
of each matrix are not uniform.
Let
maxIters = maximum number of re-estimation iterations
iters = 0
oldLogProb = .
2. The -pass
// compute
0
(i)
c
0
= 0
for i = 0 to N 1

0
(i) = (i)b
i
(O
0
)
c
0
= c
0
+
0
(i)
next i
// scale the
0
(i)
c
0
= 1/c
0
13
for i = 0 to N 1

0
(i) = c
0

0
(i)
next i
// compute
t
(i)
for t = 1 to T 1
c
t
= 0
for i = 0 to N 1

t
(i) = 0
for j = 0 to N 1

t
(i) =
t
(i) +
t1
(j)a
ji
next j

t
(i) =
t
(i)b
i
(O
t
)
c
t
= c
t
+
t
(i)
next i
// scale
t
(i)
c
t
= 1/c
t
for i = 0 to N 1

t
(i) = c
t

t
(i)
next i
next t
3. The -pass
// Let
T1
(i) = 1 scaled by c
T1
for i = 0 to N 1

T1
(i) = c
T1
next i
// -pass
for t = T 2 to 0 by 1
for i = 0 to N 1

t
(i) = 0
for j = 0 to N 1

t
(i) =
t
(i) +a
ij
b
j
(O
t+1
)
T+1
(j)
next j
// scale
t
(i) with same scale factor as
t
(i)

t
(i) = c
t

t
(i)
next i
14
next t
4. Compute
t
(i, j) and
t
(i)
for t = 0 to T 2
denom = 0
for i = 0 to N 1
for j = 0 to N 1
denom = denom +
t
(i)a
ij
b
j
(O
t+1
)
t+1
(j)
next j
next i
for i = 0 to N 1

t
(i) = 0
for j = 0 to N 1

t
(i, j) = (
t
(i)a
ij
b
j
(O
t+1
)
t+1
(j))/denom

t
(i) =
t
(i) +
t
(i, j)
next j
next i
next t
5. Re-estimate A, B and
// re-estimate
for i = 0 to N 1

i
=
0
(i)
next i
// re-estimate A
for i = 0 to N 1
for j = 0 to N 1
numer = 0
denom = 0
for t = 0 to T 2
numer = numer +
t
(i, j)
denom = denom +
t
(i)
next t
a
ij
= numer/denom
next j
next i
15
// re-estimate B
for i = 0 to N 1
for j = 0 to M 1
numer = 0
denom = 0
for t = 0 to T 2
if(O
t
== j) then
numer = numer +
t
(i)
end if
denom = denom +
t
(i)
next t
b
i
(j) = numer/denom
next j
next i
6. Compute log[P(O| )]
logProb = 0
for i = 0 to T 1
logProb = logProb + log(c
i
)
next i
logProb = logProb
7. To iterate or not to iterate, that is the question. . .
iters = iters + 1
if (iters < maxIters and logProb > oldLogProb) then
oldLogProb = logProb
goto 2
else
output = (, A, B)
end if
8 A not-so-simple example
In this section we present a classic application of Hidden Markov Models due to Cave and
Neuwirth [2]. This application nicely illustrates the strength of HMMs and has the additional
advantage that it requires no background in any specialized eld such as speech processing.
16
Suppose Marvin the Martian obtains a large body of English text, such as the Brown
Corpus [3], which totals about 1,000,000 words. Marvin, who has a working knowledge of
HMMs, but no knowledge of English, would like to determine some basic properties of this
mysterious writing system. A reasonable question he might ask is whether the characters can
be partitioned into sets so that the characters in each set are dierent in some statistically
signicant way.
Marvin might consider attempting the following. First, remove all punctuation, numbers,
etc., and convert all letters to lower case. This leaves 26 distinct letters and word space, for
a total of 27 symbols. He could then test the hypothesis that there is an underlying Markov
process (of order one) with two states. For each of these two hidden states, he assume that
the 27 symbols are observed according to xed probability distributions.
This denes an HMM with N = 2 and M = 27, where the state transition probabilities
of the A matrix and the observation probabilitiesthe B matrixare unknown, while the
the observations are the series of characters found in the text. To nd the most probable A
and B matrices, Marvin must solve Problem 3 of Section 7.
We programmed this experiment, using the rst T = 50, 000 observations (letters
converted to lower caseand word spaces) from the Brown Corpus [3]. We initialized
each element of and A randomly to approximately 1/2. The precise values used were
= [ 0.51316 0.48684 ]
and
A =
_
0.47468 0.52532
0.51656 0.48344
_
.
Each element of B was initialized to approximately 1/27. The precise values in the initial B
matrix (actually, the transpose of B) appear in the second and third columns of Figure 3.
After the initial iteration, we have
log[P(O|, )] = 165097.29
and after 100 iterations,
log[P(O| )] = 137305.28.
This shows that the model has been improved signicantly.
After 100 iterations, the model = (A, B, ) has converged to
=
_
0.00000 1.00000
_
and
A =
_
0.25596 0.74404
0.71571 0.28429
_
with the transpose of B appearing in the last two columns of Figure 3.
The most interesting part of this result is the B matrix. Without having made any
assumption about the two hidden states, the B matrix tells us that one hidden state contains
17
Initial Final
a 0.03735 0.03909 0.13845 0.00075
b 0.03408 0.03537 0.00000 0.02311
c 0.03455 0.03537 0.00062 0.05614
d 0.03828 0.03909 0.00000 0.06937
e 0.03782 0.03583 0.21404 0.00000
f 0.03922 0.03630 0.00000 0.03559
g 0.03688 0.04048 0.00081 0.02724
h 0.03408 0.03537 0.00066 0.07278
i 0.03875 0.03816 0.12275 0.00000
j 0.04062 0.03909 0.00000 0.00365
k 0.03735 0.03490 0.00182 0.00703
l 0.03968 0.03723 0.00049 0.07231
m 0.03548 0.03537 0.00000 0.03889
n 0.03735 0.03909 0.00000 0.11461
o 0.04062 0.03397 0.13156 0.00000
p 0.03595 0.03397 0.00040 0.03674
q 0.03641 0.03816 0.00000 0.00153
r 0.03408 0.03676 0.00000 0.10225
s 0.04062 0.04048 0.00000 0.11042
t 0.03548 0.03443 0.01102 0.14392
u 0.03922 0.03537 0.04508 0.00000
v 0.04062 0.03955 0.00000 0.01621
w 0.03455 0.03816 0.00000 0.02303
x 0.03595 0.03723 0.00000 0.00447
y 0.03408 0.03769 0.00019 0.02587
z 0.03408 0.03955 0.00000 0.00110
space 0.03688 0.03397 0.33211 0.01298
Figure 3: Initial and nal B transpose
the vowels while the other hidden state contains the consonants. Curiously, word-space is
more like a vowel, while y is not even sometimes a vowel. Of course, anyone familiar with
English would not be too surprised that there is a clear distinction between vowels and
consonants. But the HMM result show us that this distinction is a statistically signicant
feature inherent in the language. And, thanks to HMMs, this could easily be deduced by
Marvin the Martian, who has no other knowledge of the language.
Cave and Neuwirth [2] obtain further interesting results by allowing more than two hidden
states. In fact, they are able to obtain and sensibly interpret the results for models with up
to 12 hidden states.
18
9 Exercises
1. Write the re-estimation formulae (9), (10) and (11) directly in terms of and .
2. In the re-estimation formulae obtained in Exercise 1, substitute and

for and ,
respectively, and show that the resulting re-estimation formulae are exact.
3. Instead of using c
t
to scale the
t
(i), scale each
t
(i) by
d
t
=
1
N1

j=0

t
(j)
where the denition of

is similar to that of .
a. Using the scaling factors c
t
and d
t
show that the re-estimation formulae obtained
in Exercise 1 are exact with and

in place of and .
b. Write log[P(O| )] in terms of c
t
and d
t
.
4. Consider an HMM where for each t the state transition matrix is time dependent.
Then for each t, there is an N N row-stochastic A
t
= {a
t
ij
} that is used in place of A
in the HMM computations. For such an HMM,
a. Give pseudo-code to solve Problem 1.
b. Give pseudo-code to solve Problem 2.
5. Consider an HMM of order two, that is, an HMM where the underlying Markov process
is of order two. Then the the state at time t depends on the states at time t1 and t2
and a xed set of probabilities. For an HMM of order two,
a. Give pseudo-code to solve Problem 1.
b. Give pseudo-code to solve Problem 2.
c. Give pseudo-code to solve Problem 3.
6. Write an HMM program to solve the English language problem discussed in Section 8
with
a. Three hidden states. Explain your results.
b. Four hidden states. Explain your results.
7. Write an HMM program to solve the problem discussed in Section 8, replacing English
text with
a. French.
19
b. Russian.
c.

Chinese.
8. Perform an HMM analysis similar to that discussed in Section 8, replacing English
with Hamptonese; see
http://www.cs.sjsu.edu/faculty/stamp/Hampton/hampton.html
for information on Hamptonese.
References
[1] L. R. Rabiner, A tutorial on hidden Markov models and selected applications
in speech recognition, Proceedings of the IEEE, Vol. 77, No. 2, February 1989,
http://www.cs.ucsb.edu/~cs281b/papers/HMMs%20-%20Rabiner.pdf
[2] R. L. Cave and L. P. Neuwirth, Hidden Markov models for English, in J. D. Ferguson,
editor, Hidden Markov Models for Speech, IDA-CRD, Princeton, NJ, October 1980.
[3] The Brown Corpus of Standard American English, available for download at
http://www.cs.toronto.edu/~gpenn/csc401/a1res.html\verb
20