# 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. – Wiliam Aug 8 '18 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]]} -> #[] &,
Tally[toCount]];
toNormalize = GatherBy[toGather, #[[1, 1]] &];
SparseArray[
Flatten@Map[
With[{tot = 1/Plus @@ #[[All, 2]]},
Map[#[] -> #[] 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}} – Wiliam Aug 8 '18 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. – Hector Aug 8 '18 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)? – Wiliam Aug 8 '18 at 12:00
• I am using 8.0 ver – Wiliam Aug 8 '18 at 12:01
• @William I'll edit the answer to make it more clear. – Hector Aug 8 '18 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[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[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[[ordering]], {2}];
table = {Row@states[[{##}]], ## & @@
Table[Probability[p == s \[Conditioned] p == # && p == #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. – Wiliam Aug 8 '18 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 '18 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. – Wiliam Aug 8 '18 at 10:51
• @William, please see the update. – kglr Aug 8 '18 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 – Wiliam Aug 8 '18 at 12:37

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); // 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]


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} *)
`