Constructing higher order transition probability matrix

Question

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.

Answer 1

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} $$

Answer 2

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}$

Original answer:

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}$

Answer 3

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

Answer 4

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"]}, 
     Internal`WithLocalSettings[
      (*switch to additive assembly*)
      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]