# DiscreteMarkovProcess: preparing the transition matrix and unexpected negative results

I would like to implement a Markov process. It consists on a walk on a one-dimensional lattice, with nodes spaced $$\epsilon$$ from $$-2$$ to $$2$$. The transition probabilities for $$x \neq y$$ are as follows:

$$$$p(x,y) = \begin{cases} \frac{1}{2} \,\ \text{if} \,\, y= x\pm \epsilon \,\, \text{and} \,\,f(y) < f(x) \\ \frac{1}{2} \exp \Big(-[f(y)-f(x)] \,/ \beta\Big) \,\ \text{if} \,\, y= x\pm \epsilon \,\, \text{and} \,\,f(y) > f(x) \\ 0 \,\,\ \text{else} \end{cases}\,.$$$$

where $$\beta$$ is a given number and $$f(x) = \frac{x^4}{4}-\frac{x^2}{2}$$. For $$x=y$$, probabilites will be computed by ensuring row sums equal to $$1$$.

I constructed the probability matrix as follows:

f[x_]:= (1/4) x^4- (1/2) x^2
neighRight[beta_, eps_] :=
Table[If[f[-2 + (i)*eps] > f[-2 + (i - 1)*eps], (1/2)*
Exp[-(f[-2 + (i)*eps] - f[-2 + (i - 1)*eps])/beta],  1/2 ], {i,
2*2/eps}]
neighLeft[beta_, eps_] :=
Table[If[f[-2 + (i - 1)*eps] > f[-2 + (i)*eps], (1/2)*
Exp[-(f[-2 + (i - 1)*eps] - f[-2 + (i)*eps])/beta],  1/2 ], {i,
2*2/eps}]
diagonal[beta_, eps_] :=
Join[{1/2},
Table[1 - neighRight[beta, eps][[i + 1]] -
neighLeft[beta, eps][[i]], {i, 2*2/eps - 1}], {1/2}]
TransMatrix[beta_, eps_] :=
SparseArray[{Band[{1, 1}] -> diagonal[beta, eps],
Band[{2, 1}] -> neighLeft[beta, eps],
Band[{1, 2}] -> neighRight[beta, eps]}, {1 + 2*2/eps,
1 + 2*2/eps}];


So, the functions  neighRight and  neighLeft provide the elements to the right and left respectively of the main diagonal, which is then computed via the homonymous function using the condition that the sum of the values over the row should equal $$1$$. The transition matrix is tridiagonal so  SparseArray comes in handy.

I am sure there are better, more efficient and elegant, ways to do it and I would be grateful for any suggestions on how to assemble the transition matrix.

Nevertheless, it seems to give me the result I am looking far, after quite some "manual" checking.

Next I generate start and target distributions, such as

start = ConstantArray[0, 41] ;
end = ConstantArray[0, 41] ;
start[[31]] = 1;
end[[11]] = 1;


again, surely a caveman man to do it, but I wish the process to start on node $$31$$ and I am interested in the mean hitting time to node $$11$$.

  A = DiscreteMarkovProcess[start, TransMatrix[0.06, 0.1]];
EscTime = FirstPassageTimeDistribution[A, end];
Mean[EscTime]
3.9774*10^18


which seems reasonable. The fact, is I lower the first argument of the TransMatrix function, called $$\beta$$ above, by just a tad, e.g.

  A = DiscreteMarkovProcess[start, TransMatrix[0.06, 0.1]];
EscTime = FirstPassageTimeDistribution[A, end];
Mean[EscTime]
-4.54097*10^16


I get a negative mean hitting time, which should not be I believe. Ultimately the transition matrix seems perfectly legitimate, cannot understand what the issue is.

Could it be something to do with the very small numbers generated by the exponentials, whose argument is divided by $$\beta$$, as the latter is lowered plenty? I am at loss as to how to fix this, thanks.

• You have a missing bracket in the second row of your formula for $p(x,y)$. Could you clarify that. Jul 30, 2021 at 14:57
• @flinty, done thanks Jul 30, 2021 at 15:19
• Your code does not include the definition for f Jul 30, 2021 at 15:39
• It seems like an issue with the working precision, If you try, for example DiscreteMarkovProcess[start, SetPrecision[TransMatrix[0.05, 0.1], 20]] it gives a reasonable result. It seems like you need arbitrary precision computation for this to avoid rounding errors. Jul 30, 2021 at 15:51
• FirstPassageTimeDistribution with inexact numbers sometimes gives probabilities far outside the range 0 to 1. This has been an issue since version 10.0 and continuing with 12.3.1. High precision or exact numbers are work arounds.
Jul 30, 2021 at 17:35

I've rewritten the TransMatrix function to make it a bit easier to understand I hope:

TransMatrix[β_, ϵ_, f_, p_] := Module[{points, g, mtx},
points = Range[-2, 2, ϵ];
(* connect points up into a linear graph *)
g = RelationGraph[Abs[#1 - #2] == ϵ &, points, DirectedEdges->True];
(* set edge weights *)
g = With[{el = EdgeList[g]},
Graph[el, EdgeWeight -> (# -> p[#[[1]], #[[2]], β, ϵ, f] & /@ el)]];
(* get the weight matrix *)
(* need to make up the difference for cases when x == y if row is not normalized *)
mtx + DiagonalMatrix[1 - Total /@ mtx]
]

f[x_] := x^4/4 - x^2/2

p[x_, y_, β_, ϵ_, f_] := If[Abs[y - x] != ϵ, 0,
Piecewise[{
{1/2, f[y] < f[x]},
{1/2 Exp[(f[x] - f[y])/β], f[y] > f[x]}}, 0]]

(* use rationals here for best precision *)
trmtx = TransMatrix[6/100, 1/10, f, p];
MatrixPlot[trmtx]
(* use a high precision on trmtx *)
A = DiscreteMarkovProcess[31, SetPrecision[trmtx, 80]];
EscTime = FirstPassageTimeDistribution[A, 11];
Mean[EscTime]
(* 4491.3691... *)


This is a bit higher than experiment but that may be because it needs longer paths than 10000 steps which will take longer to compute:

SeedRandom[123];
(* mean first position of '11' in 100 random paths of max length 10000 *)
Mean[
DeleteCases[
First[FirstPosition[#, 11]] & /@ (RandomFunction[A, {0, 10000},
100]["ValueList"]), "NotFound"]
] // N

(* 3340.07 *)

• thanks a lot that is extremely instructive. Will need some time to understand exactly what you did, especially with the graph definitions. But what puzzles me the most is how different, by order of magnitudes, results are from my code. Thanks a lot. Jul 30, 2021 at 17:27
• your transition matrix is different from mine. For example, on row $3$ there are two only nonzero entries, both equal to $0.5$. This contradicts the definition. But It will take me a while to understand all your code, very instructive. Jul 30, 2021 at 18:42
• @Smerdjakov Row 3 would correspond to: {... p[-1.8, -1.9, 0.06, .1, f], p[-1.8, -1.8, 0.06, .1, f], p[-1.8, -1.7, 0.06, .1, f] ... }  giving {..., 0.5, 0, 0.5,...} according to the definition, as -1.8 is the 3rd element of Range[-2,2,.1] Jul 30, 2021 at 21:54
• I must have made a mistake but I cannot even spot it. Taking the -1.8 point, $f(-1.9) > f(-1.8) > f(-1.7)$, so the probability transition from -1.8 to -1.7 is 1/2, as per first line in the definition, but the transition from -1.8 to -1.9 will involve the exponential as per second line, and will not equal 1/2. Jul 31, 2021 at 6:23
• @Smerdjakov must be this: p[-1.9, -1.8, 0.06, .1, f] is 1/2 but p[-1.8, -1.9, 0.06, .1, f] is 0.000282953. Adding DirectedEdges -> True to the RelationGraph` fixes this. Jul 31, 2021 at 10:10