Hidden Markov Models (6)

I've written a Python script to make an HMM for the "occasionally dishonest casino." There is nothing remarkable about the code. Here is a link.

It uses the MarkovHelper module from here, modified a bit. Because long chains yield very low likelihoods, there is a function to convert the probabilities to logs:

def logs(self): def log2(x): return log(x)/log(2) fL = [log2(n) for n in self.fairProb] lL = [log2(n) for n in self.loadProb] D = dict() for k in self.tranProb: D[k] = log2(self.tranProb[k]) return fL, lL, D

There are two points where you can make edits to get more verbose output. Here is an example of output from a model with both transition probabilities set to 0.2. The top line in each group is the sequence of states, then the observed data, followed by the results from the Viterbi algorithm.

fraction correct: 0.6692 0 .. 50 LLFFFFFFFFFFFFFFLLLLFFFFFFFFFFFLLLLFFFFLFLLFLLLLLF 65466142352514466666314623515236466525163311362664 LLLLLFFFFFFFFFFLLLLLFFFFFFFFFFFLLLLFFFFFFFFFFLLLLF 50 .. 100 LLLLLLLLLFFFFFFFLLLLFFLLFFFLLLLFFFFLLFFFFLLLLLFFFF 43264622635463146366241654162162245146442626631643 FFFFFFFFFFFFFFFFLLLLFFFFFFFFFFFFFFFFFFFFFLLLLLLLLL 100 .. 150 LLLLLLFFLLLLFFFLLLFFFFFFLLLFLLLFFFLLLLFFFLLLLLLLLF 63421525664134456562136112561166465362562126526622 LFFFFFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLLLLLLL 150 .. 200 LLLLLLLLLLLLLLLLFLLLLLLFLLLLLLLFFFFFFFFFFFFLLLLFFF 66632561666263461461635121666361266546114312626611 LLLLLLLLLLLLLLLLLLLLLFFFFFLLLLLLLLLLLLFFFFFFLLLLLL

At first, this looks promising. But remember, we tell the model what all the probabilities are. So, knowing that both states are equally likely, random guessing would give 50% correct. For transition probabilities of 0.2 and 0.1, the HMM scores 0.72, 6 points higher than random. With transition probabilities of 0.2 and 0.01, the model score 0.952 correct. At least that is no worse than random!

Hidden Markov Models (5)

Building on the prrevious post in this series, we have an HMM with two states G (good) and B (bad) weather, that have characteristic emission and transition probabilities. Let's add a begin state to the model, with P = 0.5 for each possible transition.

Now, use the Viterbi algorithm to compute the most probable sequence of hidden states corresponding to some observed data. Suppose the observations are Sun-Rain-Sun. We begin with Sun. Make a table with columns for the observations and rows for the states. For the first column and row, we compute the product of the transition probability times an emission probability. This product is proportional to the probability of reaching the state, but it is not a real probability since the combined values don't add up to 1 (these are likelihoods).

Now, we add the second observation, Rain, in the second column.

Work right to left. For the G state there are two possible paths to get there, either from the previous G or from the previous B. Compute first the transition probability for B -> G times the previous cumulative probability for state B. Also compute the transition probability for G -> G times the previous cumulative probability for state G. Compare the two and take the maximum, that is, G -> G, multiply by the emission probability (Rain from G) and put the result at the bottom of the box.

Draw a horizontal arrow to remember how we got here (the traceback). Continue with the other possible state for this observation, in the second row, second column.

Notice the correspondence with the first row: probabilities on the left of the two terms in parentheses add to 1, while the cumulative probabilities being carried forward are the unchanged. In this case, the maximum value corresponds to G -> B, so we draw a diagonal arrow.

And so on:



Now, look at the right-hand column and pick the state with the highest cumulative likelihood. Follow the traceback to figure out how we got here: GBB is the most probable sequence of hidden states.

[ Two math errors were spotted by a sharp-eyed reader that alter the results in the last part, and change the most probable sequence to GBG. Sorry for the confusion. I don't have the slides anymore to edit them easily so I'll leave it as it is. ]

Hidden Markov Models (4)

In the last post I described a Python class that we will use in the future to explore the dynamic programming algorithms that are important with HMMs. But before we continue with the Python, I want to go through an example of the first of these methods, the Viterbi algorithm, which is named for Andrew Viterbi. It is applied in the case where we have settled on a model and its characteristic probabilities (the parameters), and we want to calculate:

• the most likely sequence of states given some observed output

The example I'm going to use is the weather, but with a twist from the previous example of weather by telephone. Here the weather will be the observed data. The underlying hidden states of the model are two: G (a propensity for good weather) and B (bad weather is likely). G and B can both emit either sun or rain, with characteristic emission probabilities. At each step in the chain, the states can also switch with some characteristic transition probability. (Or, stay the same, with P = 1 - P(T)).



Below is a series of observations.



For a series of n observations, in an HMM with 2 states, there are 2**n possible state sequences, each of which has some probability of emitting the observed data. Given how unbalanced the emission probabilities are for our current model, it is pretty obvious that the first state sequence above is more likely than the second one. Now, if the chain is short enough, we can calculate the probabilities for each possible sequence exactly. For example:



We observe Sun-Rain-Sun. We calculate:

GGG 0.90 * 0.60 * 0.10 * 0.60 * 0.90 = 0.02916 GGB 0.90 * 0.60 * 0.10 * 0.40 * 0.20 = 0.00432 GBG 0.90 * 0.40 * 0.80 * 0.30 * 0.90 = 0.07776 GBB 0.90 * 0.40 * 0.80 * 0.70 * 0.20 = 0.04032 BGG 0.20 * 0.30 * 0.10 * 0.60 * 0.90 = 0.00324 BGB 0.20 * 0.30 * 0.10 * 0.40 * 0.20 = 0.00048 BBG 0.20 * 0.70 * 0.80 * 0.30 * 0.90 = 0.03024 BBB 0.20 * 0.70 * 0.80 * 0.70 * 0.20 = 0.01568

So, as one would expect, GBG is most likely. In a real problem, there will be more states, and the probabilities will not be so skewed. That calls for the Viterbi algorithm. I'll show it (for this problem) in the next post.

Hidden Markov Models (3)

This post is about the "occasionally dishonest casino" from Durbin et al., Biological Sequence Analysis. To save effort, I'll just use my slide from class:



We're going to model this system in Python, starting with two sets of dice with the requisite probabilities, and extending on into HMM analysis. Here is typical output for 100 rolls of the dice. We show the state above, and the output below:

FFFFFFLFLLLFFFFFFFFFFFFLFFFFFFLLLFFFLLLFLLLLLLLLLL 44434462641561516315232623134536312512613662626653 LLLFFFFFFFFFFFLLLFFFFFFFFFFFFFFFFFFFFFFLLLLLLLLLFF 56526642512146625665112531415535611163366666636225

Following, in segments, is the code for the Model class.

import random outcomes = '123456' class Model: def __init__(self, pFL, pLF): self.emit = list(outcomes) self.fairProb = [1.0/6] * 6 self.loadProb = [0.1] * 5 + [0.5] self.tranProb = { 'F':pFL,'L':pLF } def setState(self,which=None): if which: self.state = which else: self.state = random.choice('FL')

The initialization function is kept simple. It requires a value for the transition probability to switch states in either direction. The emission probabilities are hard-coded. The state will be set at the beginning of a run. Rolling the dice (below) is straightforward. We test the result of random.random() against the cumulative probability and emit the corresponding outcome.

def rollDice(self): if self.state == 'F': L = self.fairProb else: L = self.loadProb f = random.random() S = 0 for i in range(len(L)): S += L[i] if f < S: break return self.emit[i] def transit(self): def switch(): if self.state == 'F': self.state = 'L' else: self.state = 'F' p = self.tranProb[self.state] f = random.random() if f < p: switch()

random.random() is used a second time to decide whether to switch on a particular transition. The object can generate sequences of as many states as needed, with a default of 50.

def sequence(self,N=50): rolls = list() states = list() for i in range(N): rolls.append(self.rollDice()) states.append(self.state) self.transit() return ''.join(rolls),''.join(states) if __name__ == '__main__': m = Model(pFL=0.1,pLF=0.2) m.setState('F') L = [m.sequence() for i in range(2)] for r,s in L: print s + '\n' + r + '\n' import diceStats as St St.checkEStats(m) St.checkTStats(m)

It is useful to check that model performs as expected. The code for this is in another module, which we import here. But first, to motivate us, here is typical output for this model:

1 F 0.17 L 0.097 2 F 0.166 L 0.107 3 F 0.164 L 0.099 4 F 0.162 L 0.097 5 F 0.174 L 0.096 6 F 0.164 L 0.505 F 6748 0.098 L 3252 0.204

The different outcomes 1 through 6 are emitted with the correct probabilities. The last two lines list the number of times out of 10000 that the model was in the fair or loaded state. The ratio is determined by ratio of the transition probabilities, which are also determined.

Here is the diceStats module:

outcomes = '123456' def checkEStats(m): N = 10000 m.state = 'F' fairL = [m.rollDice() for i in range(N)] m.state = 'L' loadL = [m.rollDice() for i in range(N)] for o in outcomes: print o, ' F', print str(round(fairL.count(o) * 1.0 / N, 3)).ljust(5), print ' L', print str(round(loadL.count(o) * 1.0 / N, 3)).ljust(5) print def checkTStats(m): N = 10000 L = list() for i in range(N): L.append(m.state) m.transit() totalF = L.count('F') totalL = L.count('L') D = { 'F':0,'L':0 } c1 = L[0] for c2 in L[1:]: if c1 != c2: D[c1] += 1 c1 = c2 print 'F', totalF, print str(round(D['F'] * 1.0 / totalF, 3)).ljust(5) print 'L', totalL, print str(round(D['L'] * 1.0 / totalL, 3)).ljust(5)

Hidden Markov Models (2)

A couple of weeks ago I put up an introductory post about HMMs, and promised more. I got distracted for a bit, but now I'm ready to work on it again. Before we start, let's talk a little more generally about what they are about.

Here is Andrey Markov.


Working from back to front, an HMM is

• a model
• that forms a Markov "chain"
• where the model is hidden

The Markov "property" is that we have a chain of "states" of the model and a parallel set of observations. Any particular set of observations is linear, and for each of these there is a corresponding state. The next state to be visited depends only on the current state (or a small number of recent states). And it's hidden---we don't observe the state of the model directly.

This is very cryptic, how about a concrete example?

The weather has a long history in HMM explanations (according to Dave Swofford's slides that I learned so much from, he heard about it from from some people named Rune Lyngsø (I'm guessing: this guy) via Bjarne Knudsen (bio on this page)). To get our feet wet, let's play a game called "weather by telephone."

Imagine that you and your mother live in different places, but you are a dutiful son. You call every day and she tells you about her activities. She does one of three things:

• walk
• shop at the indoor mall next door
• play scrabble.

She hates to talk about the weather, but she tells you what she's done today, and you try to guess the weather. The observed values are the chain of walk-shop-scrabble. The HMM is the weather, which we suppose has two states: sun and rain. You have to guess the state.



The model has probabilities for emission and transition. For example, if it's sunny the probability is very high that your mom will walk, whereas if it's rainy, she's most likely to stay home and play scrabble. Also, one sunny day may be more likely to be followed by another sunny day, rather than change states. We have states with transition probabilities connecting them, and for any given current state, it emits an observation drawn from a number of possible choices governed by emission probabilities that are state-specific.

Now that you're all comfortable with this analogy, I'm going to switch things around in a later post---the weather will be observed and depend on some underlying hidden state. But in the meantime, you can go back to the first post and see whether this makes sense to you in terms of the model for a functional transcription factor binding site, and the observed sequence in a particular candidate site.

Hidden Markov Models (1)

I want to start a series of posts about Hidden Markov Models or HMMs. In the spirit of the blog, these will be reports from someone who is a biologist by training, who struggled a bit with the mathematical ideas, and then found his way to a basic understanding. Nothing fancy, just an explanation that makes sense to me. And some Python code to go with it.

The best introductions that I know of are (unfortunately not free):

• The first 80 pages of Durbin
• A short article by Sean Eddy (PMID 15470472)
• A series of Dave Swofford lectures that got disappeared off the web
• A chapter by Anders Krogh in this book

Before we get into HMMs let's consider an example from Krogh explaining about Position-Specific Scoring Matrices or PSSMs. Here is a set of 5 aligned sequences:



Assume we have evidence that these are sites with the same biological function (binding data, mutations, or whatever). We want to find more examples in new sequence that we would predict have the same function.

We could use regular expressions to search for this motif. Krogh uses this example of a regular expression to match either C. elegans or Caenorhabditis elegans

"C[\.a-z]* elegans"

The notation formalizes the logic that we require first capital C, then some stuff, then a space, followed by elegans. The stuff in the middle is either a period or any lower case letter from a to z, of any length (that's the * wild-card). The backslash is just to remind ourselves that we really do mean a period.

For the sequence example, we would use [AT] [CG] [AC] [ACGT]* A [TG] [GC].

But, one problem with this approach is that we haven't used the frequency information. So, for example, our RE does not distinguish the highly implausible sequence TGCT--AGG from the consensus sequence ACAC--ATC.

Here is a PSSM for the first three positions in the sequence alignment:



We could do something similar with the last three. But a PSSM has no way simple way to deal with the variable spacing. Particularly if there are multiple variable spaces. This leads us to a simple HMM:



We proceed through a series of states which are considered to "emit" nucleotides with different frequencies, resulting in a probability for emission of any particular sequence. The states may have transition probabilities as well. In more complicated models a given sequence may be emitted by more than one series of states.

Using the basic model, we can calculate expected frequencies for different sequences. For example, the probability of the sequence ACAC--AGC is

0.8 x 0.8 x 0.8 x 0.6 x 0.4 x 0.6 x 1.0 x 0.2 x 0.8 = 0.0118

while for the second sequence TGCT--AGG we get

0.2 x 0.2 x 0.2 x 0.6 x 0.2 x 0.6 x 1.0 x 0.2 x 0.2 = 0.000023

It's probability is (satisfyingly) 500 times smaller than the first one. As you know, multiplication is tedious for us (and slow as well as error-prone for very small numbers, even on a computer). So typically the model's probabilities are specified as logarithms. And we would usually compare the probabilities according to the model with a second model based on random nucleotides at the frequencies observed in the data, although when comparing two different sequences such normalization is factored out.

You might have noticed that it was really lucky :) that both sequences had an A at the position 2 from the end. Otherwise, we'd have multiplied by zero for zero total probability. It doesn't make biological sense to do discard an otherwise excellent match for this reason. It is probably a consequence of having a limited set of authentic sequences. We say that the model is over-fitted to the data, and adjust the data using pseudocounts. So, every nucleotide that had zero probability (because it was never observed in the small set of authentic sequences) might receive a pseudo-count of 1 (or 0.1 or whatever the modeler decides is reasonable).

As an aside, PSSMs work reasonably well for many problems. I've used them a lot, and maybe I'll write about that another time.