How to solve discrete time Markov Chains in Sage in a short way

Question

I'm new to Sage.

I'm able to solve DTMC on Octave by using this short code:

a = 0.2
s = 0.6
P = [
  (1-a)*(1-a), (1-a)*a,     a*(1-a),     a*a;
  (1-a)*s,     (1-a)*(1-s), a*s,         a*(1-s);
  s*(1-a),     s*a,         (1-s)*(1-a), (1-s)*a;
  0,           s*(1-s),     (1-s)*s,     (1-s)*(1-s)+s*s;
  ]

pis = [P' - eye(size(P)); ones(1, length(P))] \ [zeros(length(P), 1); 1]

I would like to be able to do something similar in Sage. So far I'm able to solve them by using this code:

a = 0.2
s = 0.6
P = matrix(RR, 4, [
  [(1-a)*(1-a), (1-a)*a,        a*(1-a),        a*a],
  [(1-a)*s,     (1-a)*(1-s),    a*s,            a*(1-s)],
  [s*(1-a),     s*a,            (1-s)*(1-a),    (1-s)*a],
  [0,           s*(1-s),        (1-s)*s,        (1-s)*(1-s)+s*s]
  ]);
I = matrix(4, 4, 1); # I; I.parent()
s0, s1, s2, s3 = var('s0, s1, s2, s4')
eqs = vector((s0, s1, s2, s3)) * (P-I); eqs[0]; eqs[1]; eqs[2]; eqs[3]
pis = solve([
  eqs[0] == 0,
  eqs[1] == 0,
  eqs[2] == 0,
  eqs[3] == 0,
  s0+s1+s2+s3==1], s0, s1, s2, s3)

Unfortunately that code does not scale well, I have to manually edit the code to include the conditions of the equations equals to zero.

It is possible to achieve this in a way such as in Octave? It is possible to return real numbers instead of fractions?

Thanks a lot.

What do you mean by solve discrete time markov chain? Is that to find the limiting probabilities? — moldovean, May 04 '16 at 09:19

score 0 · Answer 1 · edited May 23 '17 at 12:16

I think we can create your matrices in code more analogous to what you are using in Octave, but with Python syntax. There are shortcut constructors, they just have different names.

a = 0.2
s = 0.6
P = matrix(RR, 4, [
  [(1-a)*(1-a), (1-a)*a,        a*(1-a),        a*a],
  [(1-a)*s,     (1-a)*(1-s),    a*s,            a*(1-s)],
  [s*(1-a),     s*a,            (1-s)*(1-a),    (1-s)*a],
  [0,           s*(1-s),        (1-s)*s,        (1-s)*(1-s)+s*s]
  ]);
M = (P.conjugate_transpose() - identity_matrix(P.ncols())).stack(matrix(1,P.ncols(),[1]*P.ncols()))
V = vector( [0]*P.ncols()+[1])

However, maybe I missed something in the translation; I am not very familiar with Octave, and just looked through the documentation. But the prime seems to be conjugate transpose, eye seems to be identity, and ones is all-ones, and zeros is zeros.

sage: M \ V
ValueError: matrix equation has no solutions

I note your answer here so I assume you can take it from here. Maybe there was a transpose difference in the Sage implementation of the backslash operator?

As to the fractions, this is because Maxima provides Sage's solving mechanisms, and it is more symbolically oriented - so uses fractions whenever possible. We try to keepfloat:true there as much as we can, but it doesn't always like this.

How to solve discrete time Markov Chains in Sage in a short way

1 Answers1