DiscreteMarkovProcess: preparing the transition matrix and unexpected negative results

Question

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:

\begin{equation} 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}\,. \end{equation}

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.

Answer 1

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 *)
  mtx = WeightedAdjacencyMatrix[g];
  (* 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 *)