# Problem with a Markov's Matrix

I recently discovered the Markov process.
To simplify my work, I use the following code taken from the DiscreteMarkovProcess documentation:

K = DiscreteMarkovProcess[{1, 0, 0}, ({
{0, 1/2, 1/2},
{1/2, 0, 1/2},
{1/2, 1/2, 0}
})
];

data = RandomFunction[K, {0, 10}]


However, I am facing a case where the code does not seem to work. In fact, when I use the following matrix RandomFunction does not return a result:

    {{0.666667, 0.333333, 0., 0., 0., 0., 0., 0., 0.}, {0., 0., 0., 1.,
0., 0., 0., 0., 0.}, {0.0263158, 0., 0.315789, 0.105263, 0.526316,
0.0263158, 0., 0., 0.}, {0., 0., 0., 0., 0.655172, 0.344828, 0., 0.,
0.}, {0., 0., 0.0481928, 0.060241, 0.678715, 0.176707, 0., 0.,
0.0361446}, {0., 0., 0.210526, 0.105263, 0.684211, 0., 0., 0.,
0.}, {0., 0., 0., 1., 0., 0., 0., 0., 0.}, {0., 0., 1., 0., 0., 0.,
0., 0., 0.}, {0., 0., 0., 0., 0.181818, 0.181818, 0.272727,
0.181818, 0.181818}}


To investigate the problem I wrote the following code:

    TheV = RandomReal[{0, 1}, {4, 4}];
no = Total /@ TheV;
TheVno = Table[TheV[[i]]/no[[i]], {i, 1, Length@no, 1}];
K = DiscreteMarkovProcess[1, TheVno];
data = RandomFunctionK, {0, 10}]


This code seems to show that in some cases RandomFunction does not work. What is happening in these cases ? Do I need to replace the command RandomFunction with something else ?

Edit n°1

To better investigate the problem I wrote the following code :

MarkovTest[x_, y_] :=
Block[
{step1, step2, step31, step311, step32, step322, step4,
step41, step5, step51, step6, step61, res},

step1 = Table[
RandomReal[{0, 1}, {BlockRandom[RandomInteger[{1, 18}]],
RandomInteger[{1, 18}]}],
{i, 1, x, 1}
];
step2 = Table[
Total /@ step1[[i]],
{i, 1, x, 1}
];
step31 = Table[
N[step1[[i, j]]/step2[[i, j]], y],
{i, 1, Length@step1,1}, {j, 1, Length@step1[[i]],1}
];
step311 = Table[
Rationalize@(step1[[i, j]]/step2[[i, j]]),
{i, 1, Length@step1,1}, {j, 1, Length@step1[[i]],1}
];
step32 = Table[
DiscreteMarkovProcess[1, step31[[i]]],
{i, 1, x, 1}
];
step322 = Table[
DiscreteMarkovProcess[1, step311[[i]]],
{i, 1, x, 1}
];
step4 = Table[
RandomFunction[step32[[i]], {0, 5}],
{i, 1, x, 1}
];
step41 = Table[
RandomFunction[step322[[i]], {0, 5}],
{i, 1, x, 1}
];

step6 = Count[step5, TemporalData];
step61 = Count[step51, TemporalData];

res = {step6, step61}
];


With this function, you can see the different cases where RandomFunction does not work. As mentioned by ubpdqn the main problem seems to be the arithmetic precision.
I continue to work on the problem.

PS : Obviously, I checked that the given matrix is stochastic....

Apart from rounding errors, which Mathematica normalizes, this does work.

mat = {{0.666667, 0.333333, 0., 0., 0., 0., 0., 0., 0.}, {0., 0., 0.,
1., 0., 0., 0., 0., 0.}, {0.0263158, 0., 0.315789, 0.105263,
0.526316, 0.0263158, 0., 0., 0.}, {0., 0., 0., 0., 0.655172,
0.344828, 0., 0., 0.}, {0., 0., 0.0481928, 0.060241, 0.678715,
0.176707, 0., 0., 0.0361446}, {0., 0., 0.210526, 0.105263,
0.684211, 0., 0., 0., 0.}, {0., 0., 0., 1., 0., 0., 0., 0.,
0.}, {0., 0., 1., 0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0.,
0.181818, 0.181818, 0.272727, 0.181818, 0.181818}};
dmk = Quiet@DiscreteMarkovProcess[{1, 0, 0, 0, 0, 0, 0, 0, 0}, mat];
Graph[dmk]


I leave it up to OP to check whether this is intended process.

Using RandomFunction:

llp = ListLinePlot[#,
InterpolationOrder -> 0] & /@ (RandomFunction[dmk, {0, 20}, 100][
"Paths"]);


The following gif was generated from llp:

• I was rather surprised that in your case the command RandomFunction works with my matrix. However, this confirms that the problem seems to come from the management of the arithmetic precision by Mathematica. – Doedalos Nov 10 '14 at 14:18

Your random matrix is not a valid transition matrix:

mat = {{0.666667, 0.333333, 0., 0., 0., 0., 0., 0., 0.}, {0., 0., 0., 1., 0., 0., 0., 0., 0.}, {0.0263158, 0., 0.315789, 0.105263,
0.526316, 0.0263158, 0., 0., 0.}, {0., 0., 0., 0., 0.655172, 0.344828, 0., 0., 0.}, {0., 0., 0.0481928, 0.060241, 0.678715,
0.176707, 0., 0., 0.0361446}, {0., 0., 0.210526, 0.105263, 0.684211, 0., 0., 0., 0.}, {0., 0., 0., 1., 0., 0., 0., 0., 0.}, {0.,
0., 1., 0., 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0.181818, 0.181818, 0.272727, 0.181818, 0.181818}};

Total[mat[[All, 1 ;; -1]]]


You can see the probabilities do not add to 1.

While the first matrix you used which work does:

mat = {{0, 1/2, 1/2}, {1/2, 0, 1/2}, {1/2, 1/2, 0}};
Total[mat[[All, 1 ;; -1]]]
(* {1, 1, 1} *)


You can't just use any random matrix for DiscreteMarkovProcess, the matrix must have certain properties.

http://en.wikipedia.org/wiki/Stochastic_matrix

• The rows should add up to 1, not the columns. – Sjoerd C. de Vries Nov 9 '14 at 16:37
• Which means that your answer is incorrect. I suggest you delete it. I don't understand the upvotes here. – Sjoerd C. de Vries Nov 11 '14 at 10:47

After some research, I found a reliable way to always get the results from the command RandomFunction.

A stochastic matrix is a real square matrix, with each row summing to 1.
Even if this condition is satisfied, it often happens that the command RandomFunction does not work. This is due to how Mathematica manages the arithmetic accuracy of these numbers.

You can see a basic solution through the following code :

    l1 = RandomReal[{0, 1}, {BlockRandom[RandomInteger[{10, 15}]],
RandomInteger[{10, 15}]}];
l2 = Total /@ l1;
l3 = Table[Round[l1[[i]]/l2[[i]], 0.0005], {i, 1, Length@l1, 1}];
l4 = Map[Rationalize[#] &, l3, {2}];

Head@(RandomFunction[DiscreteMarkovProcess[1, l4], {0, 5}]) == TemporalData


In the previous lines of code, even if Mathematica must make a correction, it still provide TemporalData.

To avoid possible problems with matrix in the Markov process, we can use the following function :

MarkovGarantee[x_] := Map[Rationalize[Round[#, 0.00005]] &, x, {2}]


This solution is empirical. I'm still open to a more elegant solutions.

Edit n°1

A better solution :

MarkovGarantee[x_] :=
Block[
{step1, step2, res},

step1 = Map[Rationalize[Round[#, 0.00005]] &, x, {2}];
step2 = (1 - Total /@ step1)/Length@x;
res = Table[step1[[i]] + step2[[i]], {i, 1, Length@step1, 1}]
];


This solution don't generate any messages of error.

Edit n°2

A better solution :

MarkovGarantee2[x_] :=
Block[
{step1, step2, step3, step4, step5, step6, step7},

step1 = Map[Rationalize[Round[#, 0.00005]] &, x, {2}];
step2 = (1 - Total /@ step1)/Length@x;
step3 = step2*(Length@x);
step4 = Map[Count[#, Except[0]] &, step1];
step5 = Map[step3[[#]]/step4[[#]] &, Range[1, Length@step3, 1]];
step6 = Position[#, y_ /; y != 0] & /@ step1;
step7 =
Table[
MapAt[# + step5[[i]] &, step1[[i]], step6[[i]]],
{i, 1, Length@step6, 1}
]
];


This solution enable to avoid the appearance of negative number.