Finding the likeliest path in a Markov process

Question

With Mathematica 9, we have the addition of various processes, among which the discrete Markov process. Given a transition probability matrix m, such a process is defined as follows:

m = {{1/4, 1/4, 1/4, 1/4, 0, 0},
     {1/3, 0, 1/3, 1/3, 0, 0},
     {0, 1/2, 0, 1/4, 1/16, 3/16},
     {0, 0, 0, 1/2, 1/2, 0},
     {0, 0, 0, 0, 1/2, 1/2},
     {0, 0, 0, 0, 0, 1}
    };

proc = DiscreteMarkovProcess[1, m];

g1 = Graph[proc]

Various properties of a given Markov process can be found using the function MarkovProcessProperties. Probability can be used (among other things) to answer questions about the likelihood to end up in a certain state after a given number of steps.

However, I did not find a function that yields the path that is most likely to be taken between a given a starting and ending state. How would one find this path?

Answer 1

This sounds very similar to a traveling salesman problem for which we have the built-in (as of v8) FindShortestPath function. However, this function minimizes the overall path length (which is the sum of all sub-paths), whereas we need a functions that maximizes the path probability (which is the product of the transition probabilities along the path).

Since maximizing a product is equivalent to minimizing the sum of the negated logarithms this is easy enough. In this case:

g2 = WeightedAdjacencyGraph[-Log[m]]

sp2= FindShortestPath[g2, 1, 6]

{1, 4, 5, 6}

HighlightGraph[g1, sp2]

The shortest path in terms of the original matrix weights would have been:

sp1 = FindShortestPath[g1, 1, 6]

{1, 3, 6}.

Test that the former is indeed more likely than the latter:

pathProbability[m_?MatrixQ, path_List] := Times @@ (m[[##]] & @@@ Partition[path, 2, 1])

pathProbability[m, #] & /@ {sp1, sp2}

{3/64, 1/16}

---

UPDATE

Since version 10 Mathematica has acquired tools for hidden Markov processes. The likeliest path can now be calculated by setting the starting and ending state to 100% emitting and all other states to hidden, non-emitting states. FindHiddenMarkovStates will now do the job for us using Viterbi decoding, though it doesn't really yield shorter or more elegant code. It's also a lot (by a factor of more than 10) slower.

hmm = HiddenMarkovProcess[1, m, {{1, 0, 0, 0, 0, 0}, None, None, None, None, {0, 0, 0, 0, 0, 1}}];
obs = {1, 6};
FindHiddenMarkovStates[obs, hmm]

{1, 4, 5, 6}