35

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]

Mathematica graphics

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?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Sjoerd C. de Vries
  • 65,815
  • 14
  • 188
  • 323
  • 2
    Isn't that just what Viterbi's algorithm does? – Niki Estner Feb 12 '13 at 08:22
  • @nikie I wasn't aware of that algorithm. Thanks. According to Wikipedia this is used for hidden markov models, but I feel that that implies it should be easily convertible to MM as well. The implementation will take a few more lines of code than the one given below, which I believe is already optimal. I may be wrong here, though. – Sjoerd C. de Vries Feb 12 '13 at 09:26

1 Answers1

34

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]]

Mathematica graphics

sp2= FindShortestPath[g2, 1, 6]

{1, 4, 5, 6}

HighlightGraph[g1, sp2]

Mathematica graphics

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}

Sjoerd C. de Vries
  • 65,815
  • 14
  • 188
  • 323