# Constructing higher order transition probability matrix

Recently I asked a question here about how to construct a transition probability matrix given the following list:

x = {"A", "A", "A", "E", "D", "D", "D", "C", "B", "E", "E", "E", "D",
"B", "A", "D", "B", "E", "C", "A", "D", "A", "A", "A", "A", "C",
"C", "C", "D", "D", "E"}


For which one can get the following matrix: (see the detail from the previous question)

$$\begin{array}{cccccc} & \text{A} & \text{B} & \text{C} & \text{D} & \text{E} \\ \text{A} & \frac{5}{9} & 0 & \frac{1}{9} & \frac{2}{9} & \frac{1}{9} \\ \text{B} & \frac{1}{3} & 0 & 0 & 0 & \frac{2}{3} \\ \text{C} & \frac{1}{5} & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} & 0 \\ \text{D} & \frac{1}{8} & \frac{1}{4} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} \\ \text{E} & 0 & 0 & \frac{1}{5} & \frac{2}{5} & \frac{2}{5} \\ \end{array}$$

above is equivalent of partitioning list $x$ into sublists with size 2 and offset of 1, then counting each element and divide it by the sum of the row. The command to find the right partition is Partition[x, 2, 1] (again I refer you to the previous question). Now what if we want to find the higher order transition matrix? For example the second order would be Partition[x, 3, 1] and the expected matrix shall look like:

$$\begin{array}{cccccc} & A & B & C & D &E \\ AA & P_{AA,A} & P_{AA,B} & P_{AA,C} & P_{AA,D} &P_{AA,E}\\ AB & P_{AB,A} & P_{AB,B} & P_{AB,C} & P_{AB,D} &P_{AB,E}\\ AC & P_{AC,A} & P_{AC,B} & P_{AC,C} & P_{AC,D} &P_{AC,E}\\ AD & P_{AD,A} & P_{AD,B} & P_{AD,C} & P_{AD,D} &P_{AD,E}\\ AE & P_{AE,A} & P_{AE,B} & P_{AE,C} & P_{AE,D} &P_{AE,E}\\ \vdots & \vdots & \vdots & \vdots &\vdots & \vdots\\ EC & P_{EC,A} & P_{EC,B} & P_{EC,C} & P_{EC,D} &P_{EC,E}\\ ED &P_{ED,A} & P_{ED,B} & P_{ED,C} & P_{ED,D} &P_{ED,E}\\ EE & P_{EE,A} & P_{EE,B} & P_{EE,C} & P_{EE,D} &P_{EE,E}\\ \end{array}$$

In general the dimension of the matrix follows $\{|S|^n,|S|\}$, where n is the order of the Markov chain.

• Thanks for the edit @rhermans it was really bugging me. Aug 8, 2018 at 9:59

The following code is just brute-force. But at least yields the expected results. Also, it can be used for any order.

The first parameter is the data. The second parameter is the order.

probM[data_, ord_] :=
Module[{uniques = Union[data], acc = 0, len, trans, trPre, tData,
toCount, toGather, toNormalize},
trans = Dispatch@Thread[uniques -> Range[len = Length[uniques]]];
trPre = Dispatch@Flatten[Array[{##} -> ++acc &, ConstantArray[len, ord]]];
tData = Replace[data, trans, {1}];
toCount = Partition[tData, ord + 1, 1];
toGather = Map[{Replace[#[[1, ;; -2]], trPre], #[[1, -1]]} -> #[[2]] &,
Tally[toCount]];
toNormalize = GatherBy[toGather, #[[1, 1]] &];
SparseArray[
Flatten@Map[
With[{tot = 1/Plus @@ #[[All, 2]]},
Map[#[[1]] -> #[[2]] tot &, #]] &, toNormalize]]];


Let us check the dimensions of the first three orders.

Table[probM[x, i] // Dimensions, {i, 3}]
(*{{5, 5}, {25, 5}, {125, 5}}*)


As for the efficiency of probM, I tried replacing some of the Map with ParallelMap but it did not yield any improvement. You might want to combine with niceties from the other answer. For example, use ArrayComponents instead of dispatch tables.

In any case, check the second order table:

$$\begin{array}{cccccc} \text{} & \text{A} & \text{B} & \text{C} & \text{D} & \text{E} \\ \text{AA} & \frac{3}{5} & 0 & \frac{1}{5} & 0 & \frac{1}{5} \\ \text{AB} & 0 & 0 & 0 & 0 & 0 \\ \text{AC} & 0 & 0 & 1 & 0 & 0 \\ \text{AD} & \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 \\ \text{AE} & 0 & 0 & 0 & 1 & 0 \\ \text{BA} & 0 & 0 & 0 & 1 & 0 \\ \text{BB} & 0 & 0 & 0 & 0 & 0 \\ \text{BC} & 0 & 0 & 0 & 0 & 0 \\ \text{BD} & 0 & 0 & 0 & 0 & 0 \\ \text{BE} & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} \\ \text{CA} & 0 & 0 & 0 & 1 & 0 \\ \text{CB} & 0 & 0 & 0 & 0 & 1 \\ \text{CC} & 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ \text{CD} & 0 & 0 & 0 & 1 & 0 \\ \text{CE} & 0 & 0 & 0 & 0 & 0 \\ \text{DA} & 1 & 0 & 0 & 0 & 0 \\ \text{DB} & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\ \text{DC} & 0 & 1 & 0 & 0 & 0 \\ \text{DD} & 0 & 0 & \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\ \text{DE} & 0 & 0 & 0 & 0 & 0 \\ \text{EA} & 0 & 0 & 0 & 0 & 0 \\ \text{EB} & 0 & 0 & 0 & 0 & 0 \\ \text{EC} & 1 & 0 & 0 & 0 & 0 \\ \text{ED} & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ \text{EE} & 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ \end{array}$$

• Thanks a lot, what command did you use to get the final table? when I run your code I get {{5, 5}, {22, 5}, {109, 5}} Aug 8, 2018 at 11:53
• @William Try probM[x, 2] // MatrixForm. As for your result, it seems that you run Table[probM[x, i] // Dimensions, {i, 3}]. But that should have returned {{5, 5}, {25, 5}, {125, 5}}. What version of MMA are you running? Mine is 9.0. Aug 8, 2018 at 11:57
• Perfect, and how can I change the orders (this is second order, if I want to check the code with 1, 3 and other orders)? and the offset (possibly)? Aug 8, 2018 at 12:00
• I am using 8.0 ver Aug 8, 2018 at 12:01
• @William I'll edit the answer to make it more clear. Aug 8, 2018 at 12:02

Update: Using EmpiricalDistribution and MarginalDistribution to compute the conditional probabilities:

ClearAll[transitionProb]
transitionProb[step_: 1][x_] := Module[{states = DeleteDuplicates@x,
ed = EmpiricalDistribution[Partition[ArrayComponents @ x, step + 1, 1]],
ordering, tuples, md, condpdF},
ordering = Ordering[states]; tuples = Tuples[ordering, step];
md = MarginalDistribution[ed, Range[step]];
condpdF[u__, w_] := If[PDF[md, {u}] === 0, 0, PDF[ed, {u, w}]/PDF[md, {u}]];
Prepend[{Row @ states[[{##}]],
## & @@ Table[## & @@ condpdF[##, i], {i, ordering}]} & @@@ tuples,
Prepend[states[[ordering]], ""]]]


Examples:

transitionProb[2][x] // Grid[#, Dividers -> All] & // TeXForm


$\begin{array}{|c|c|c|c|c|c|} \hline \text{} & \text{A} & \text{B} & \text{C} & \text{D} & \text{E} \\ \hline \text{A}\text{A} & \frac{3}{5} & 0 & \frac{1}{5} & 0 & \frac{1}{5} \\ \hline \text{A}\text{B} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{A}\text{C} & 0 & 0 & 1 & 0 & 0 \\ \hline \text{A}\text{D} & \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 \\ \hline \text{A}\text{E} & 0 & 0 & 0 & 1 & 0 \\ \hline \text{B}\text{A} & 0 & 0 & 0 & 1 & 0 \\ \hline \text{B}\text{B} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{B}\text{C} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{B}\text{D} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{B}\text{E} & 0 & 0 & \frac{1}{2} & 0 & \frac{1}{2} \\ \hline \text{C}\text{A} & 0 & 0 & 0 & 1 & 0 \\ \hline \text{C}\text{B} & 0 & 0 & 0 & 0 & 1 \\ \hline \text{C}\text{C} & 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ \hline \text{C}\text{D} & 0 & 0 & 0 & 1 & 0 \\ \hline \text{C}\text{E} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{D}\text{A} & 1 & 0 & 0 & 0 & 0 \\ \hline \text{D}\text{B} & \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} \\ \hline \text{D}\text{C} & 0 & 1 & 0 & 0 & 0 \\ \hline \text{D}\text{D} & 0 & 0 & \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\ \hline \text{D}\text{E} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{E}\text{A} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{E}\text{B} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{E}\text{C} & 1 & 0 & 0 & 0 & 0 \\ \hline \text{E}\text{D} & 0 & \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ \hline \text{E}\text{E} & 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ \hline \end{array}$

transitionProb[1][x] // Grid[#, Dividers -> All] & // TeXForm


$\begin{array}{|c|c|c|c|c|c|} \hline \text{} & \text{A} & \text{B} & \text{C} & \text{D} & \text{E} \\ \hline \text{A} & \frac{5}{9} & 0 & \frac{1}{9} & \frac{2}{9} & \frac{1}{9} \\ \hline \text{B} & \frac{1}{3} & 0 & 0 & 0 & \frac{2}{3} \\ \hline \text{C} & \frac{1}{5} & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} & 0 \\ \hline \text{D} & \frac{1}{8} & \frac{1}{4} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} \\ \hline \text{E} & 0 & 0 & \frac{1}{5} & \frac{2}{5} & \frac{2}{5} \\ \hline \end{array}$

states = DeleteDuplicates[x];
ordering = Ordering[states];
data = ArrayComponents@x ;
estproc = EstimatedProcess[data, DiscreteMarkovProcess[Length@states]];
tuples = Tuples[Range[5][[ordering]], {2}];
table = {Row@states[[{##}]], ## & @@
Table[Probability[p[3] == s \[Conditioned] p[1] == # && p[2] == #2,
p \[Distributed] estproc], {s, Range[Length @ states]}]} & @@@ tuples ;

TeXForm @ Grid[Prepend[table, Prepend[states[[ordering]], ""]], Dividers -> All]


$\begin{array}{|c|c|c|c|c|c|} \hline \text{} & \text{A} & \text{B} & \text{C} & \text{D} & \text{E} \\ \hline \text{AA} & \frac{5}{9} & \frac{1}{9} & \frac{2}{9} & \frac{1}{9} & 0 \\ \hline \text{AB} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{AC} & \frac{1}{5} & 0 & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} \\ \hline \text{AD} & \frac{1}{8} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} & \frac{1}{4} \\ \hline \text{AE} & 0 & \frac{2}{5} & \frac{2}{5} & \frac{1}{5} & 0 \\ \hline \text{BA} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{BB} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{BC} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{BD} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{BE} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{CA} & \frac{5}{9} & \frac{1}{9} & \frac{2}{9} & \frac{1}{9} & 0 \\ \hline \text{CB} & \frac{1}{3} & \frac{2}{3} & 0 & 0 & 0 \\ \hline \text{CC} & \frac{1}{5} & 0 & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} \\ \hline \text{CD} & \frac{1}{8} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} & \frac{1}{4} \\ \hline \text{CE} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{DA} & \frac{5}{9} & \frac{1}{9} & \frac{2}{9} & \frac{1}{9} & 0 \\ \hline \text{DB} & \frac{1}{3} & \frac{2}{3} & 0 & 0 & 0 \\ \hline \text{DC} & \frac{1}{5} & 0 & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} \\ \hline \text{DD} & \frac{1}{8} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} & \frac{1}{4} \\ \hline \text{DE} & 0 & \frac{2}{5} & \frac{2}{5} & \frac{1}{5} & 0 \\ \hline \text{EA} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{EB} & 0 & 0 & 0 & 0 & 0 \\ \hline \text{EC} & \frac{1}{5} & 0 & \frac{1}{5} & \frac{2}{5} & \frac{1}{5} \\ \hline \text{ED} & \frac{1}{8} & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} & \frac{1}{4} \\ \hline \text{EE} & 0 & \frac{2}{5} & \frac{2}{5} & \frac{1}{5} & 0 \\ \hline \end{array}$

• Thanks a lot for this, however as a check I tried 'Partition[x, 3, 1] // Counts' which shows that I have 3 AAA instead of 5, this happens all over the table my worry is that the state ordering and data which you defined are not matched, for example in ordering E has value of 5, while in data it is 2. Aug 8, 2018 at 10:28
• @William, re ordering of states, that's why we sort them using ordering so that 2 corresponds to E. Re the discrepancy between Partition[x, 3, 1] // Counts and the Prob[A|AA] in the table above, i think it is because table is based on the TransitionMatrix of estproc and estproc is based on one-step transitions (Partition[x,2,1]).
– kglr
Aug 8, 2018 at 10:39
• So basically DiscreteMarkovProcess doesn't have a memory, because in second order the transition is remembering the two steps behind that's why it is AA,A and not for example A,A,A, interesting. Aug 8, 2018 at 10:51
• @William, please see the update.
– kglr
Aug 8, 2018 at 12:26
• thank you @kglr, it indeed works fine and is shorter yet if you change the variable in Partition[data, 3, 1] for example to Partition[data, 2, 1] or else, the final table will not come out appropriate. I think hector's answer is good because you can change the order Aug 8, 2018 at 12:37

You can use CrossTensorate from the package CrossTabulate.m, which I used and referenced in my answer of the previous question.

The making of contingency tensors with that function is discussed in this blog post: "Contingency tables creation examples".

In general, though, I would say it is better to use Tries with Frequencies or nested associations.

tmat3 = CrossTensorate[Count == 1 + 2 + 3, Partition[x, 3, 1]];

tmat4 = CrossTensorate[Count == 1 + 2 + 3 + 4, Partition[x, 4, 1]];

tmat3["XTABTensor"] = #/(Total[#, {Length[Dimensions[#]]}] /. {0 -> 1}) &@tmat3["XTABTensor"];
tmat4["XTABTensor"] = #/(Total[#, {Length[Dimensions[#]]}] /. {0 -> 1}) &@tmat4["XTABTensor"];

Grid[{{"tmat3", "tmat4"}, {MatrixForm[tmat3], MatrixForm[tmat4]}}]


ArrayRules[tmat3["XTABTensor"]]

(* {{1, 1, 1} -> 3/5, {1, 1, 5} -> 1/5, {1, 5, 4} ->
1, {1, 4, 2} -> 1/2, {1, 4, 1} -> 1/2, {1, 1, 3} -> 1/
5, {1, 3, 3} -> 1, {2, 5, 5} -> 1/2, {2, 1, 4} -> 1, {2, 5, 3} -> 1/
2, {3, 2, 5} -> 1, {3, 1, 4} -> 1, {3, 3, 3} -> 1/2, {3, 3, 4} -> 1/
2, {3, 4, 4} -> 1, {4, 4, 4} -> 1/3, {4, 4, 3} -> 1/3, {4, 3, 2} ->
1, {4, 2, 1} -> 1/2, {4, 2, 5} -> 1/2, {4, 1, 1} -> 1, {4, 4, 5} ->
1/3, {5, 4, 4} -> 1/2, {5, 5, 5} -> 1/2, {5, 5, 4} -> 1/
2, {5, 4, 2} -> 1/2, {5, 3, 1} -> 1, {_, _, _} -> 0} *)


As a variant of my answer to the linked question, the following should work correctly and efficiently.

Some random data to work with:

x = RandomChoice[Alphabet["English", "IndexCharacters"], 1000000];


Creating the probability tensor P:

n = 2;
data = Flatten[ToCharacterCode[x]] - (ToCharacterCode["A"][[1]] - 1); // AbsoluteTiming // First
A = With[{spopt = SystemOptions["SparseArrayOptions"]},
InternalWithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],

(*assemble matrix*)
SparseArray[Partition[data, n + 1, 1] -> 1, ConstantArray[Max[data], n + 1] ],

(*reset "SparseArrayOptions" to previous value*)
SetSystemOptions[spopt]]]; // AbsoluteTiming // First
P = #/N[Total[Abs[#], {n + 1}] /. 0 -> 1] &@Flatten[A, n - 1];


0.717521

0.184357

The row labels of P should be

Tuples[Sort[DeleteDuplicates[x]], n]
`