1
$\begingroup$

this question is based on a previous post of mine. I am simulating a one-dimensional Markov process on a linear graph. Transition probabilities are different than zero only for nearest neighbours of the current state as well as for transition back to the current state itself. Transition probabilities depend on a parameter $\beta$, an energy function $f$ and are computed via a function $p$. Transition probabilities for a jump to the current state are computed by imposing that the sum over rows equals $1$. Following the excellent answer by Flinty, the transition matrix is constructed with

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

Setting $p$ and $f$ as

f[x_] := ((1/2)*x^2 - 1)*(1/2)*x^2;
p[x_, y_, \[Beta]_, \[Epsilon]_, f_] := 
 If[Abs[y - x] != \[Epsilon], 0, 
  Piecewise[{{1/3000, 
     f[y] < f[x]}, {1/3000 Exp[(f[x] - f[y])/\[Beta]], f[y] > f[x]}}, 
   0]];

the following transition matrix van be computed computed, taking for example $\beta = 0.25$ and $\epsilon = 1$

trmtx = TransMatrixFlint[0.25, 1, f, p]

whose sums over rows are not that close to $1$

Total[trmtx]
(*{0.999789, 1.00054, 1., 1., 1., 1., 0.999667}*)

I have tried to construct the same transition probability matrix without using the admittedly elegant graph machinery used by Flinty in his past answer, and get different results

test[\[Beta]_, \[Epsilon]_, f_, p_] := 
  Table[p[i, j, \[Beta], \[Epsilon], 
    f], {i, {0, 1, 2, 3, 4, 5, 6}}, {j, {0, 1, 2, 3, 4, 5, 6}}];
c = test[0.25`10, 1, f, p]
a = DiagonalMatrix[1 - Total[test[0.25, 1, f, p]]]
Total[c + a]
(*{1.0000000000000, 1., 1., 1., 1., 1., 1.} *)

I am frankly unable to point the reasons for this behaviour, any hint would be most appreciated, thanks.

If it could be of any convenience, I add below the tow transition matrices. Why are some figures rationalised in the second approach, and not the first? Why are some other so different?

trmtx // MatrixForm
c + a // MatrixForm
(* {{0.999667, 0.000333333, 0., 0., 0., 0., 0.}, {0.000122626, 0.999877, 
  4.11366*10^-8, 0., 0., 0., 0.}, {0., 0.000333333, 0.999667, 
  4.33194*10^-28, 0., 0., 0.}, {0., 0., 0.000333333, 0.999667, 
  3.99454*10^-74, 0., 0.}, {0., 0., 0., 0.000333333, 0.999667, 
  1.21763*10^-156, 0.}, {0., 0., 0., 0., 0.000333333, 0.999667, 
  4.63191*10^-286}, {0., 0., 0., 0., 0., 0.000333333, 0.999667}} *)
(* {{0.9998773735196, 1/3000, 0, 0, 0, 0, 0}, {0.0001226264804, 
  0.9993333333333, 4.11366014*10^-8, 0, 0, 0, 0}, {0, 1/3000, 
  0.9996666255301, 4.3319381*10^-28, 0, 0, 0}, {0, 0, 1/3000, 
  0.9996666666667, 3.9945435*10^-74, 0, 0}, {0, 0, 0, 1/3000, 
  0.9996666666667, 1.2176307*10^-156, 0}, {0, 0, 0, 0, 1/3000, 
  0.9996666666667, 4.631909*10^-286}, {0, 0, 0, 0, 0, 1/3000, 
  1.0000000000000}}
*)

EDIT

In general I always have accuracy issues with Markov Processes. For example defining a function that computes the transition matrix as above,

TransMatrixMM[\[Beta]_, \[Epsilon]_, f_, p_] = 
  test[0.25, 1, f, p] + DiagonalMatrix[1 - Total[test[0.25, 1, f, p]]];
a = DiagonalMatrix[1 - Total[test[0.25`10, 1, f, p]]];

it seems its rows nicely sum to $1$

Total[TransMatrixMM[0.25, 1., f, p]]
(*{1., 1., 1., 1., 1., 1., 1.} *)

yet, if I try to compute the mean time to state $6$ for example, starting from $1$, with

arrheniusMM[func_, temp_, eps_] := 
 Mean[FirstPassageTimeDistribution[
   DiscreteMarkovProcess[1, 
    SetPrecision[TransMatrixMM[temp, eps, func, p], 120]], 6]];
{arrheniusMM[f, 0.25, 1], arrhenius[f, 0.25, 1]} 

I get the warning

DiscreteMarkovProcess::dmnorm: The transition matrix had some row sums that were not 1, so those rows have been adjusted by dividing by the row sum or by setting the diagonal element to 1.

cannot understand why at all.

$\endgroup$
1
  • $\begingroup$ I think I understood why the two methods give different results. I missed that Total [list] sums over columns, while Total /@ list does sum over rows, as desired. I have got that bit to work, but I still get those warning on the rows not being $1$. Continue to look at it $\endgroup$
    – Smerdjakov
    Feb 5 at 19:05

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy