Add to mpmath some useful routines specialized for stochastic (or Markov) matrices
In [1]:
import mpmath as mp
In [2]:
from eigen_markov import stoch_eig, gth_solve
stoch_eig returns the stochastic eigenvector (or stationary probability distribution vector) of
an irreducible
stochastic martix P,
i.e., the solution to $x P = x$, normalized so that its 1-norm equals one.
In [3]:
P = mp.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
In [4]:
print P
In [5]:
x = mp.mp.stoch_eig(P)
In [6]:
print x
In [7]:
print x * P
In [8]:
x*P == x
Out[8]:
In [9]:
type(x[0])
Out[9]:
Fast low-precision arithmetic:
In [10]:
P_fp = mp.fp.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
x_fp = mp.fp.stoch_eig(P)
In [11]:
print x_fp
In [12]:
type(x_fp[0])
Out[12]:
Arbitrary-precision interval arithmetic:
In [13]:
P_iv = mp.iv.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
x_iv = mp.iv.stoch_eig(P_iv)
In [14]:
print x_iv
In [15]:
type(x_iv[0])
Out[15]:
In [16]:
print x_iv * P_iv
In [17]:
for value, interval in zip(x, x_iv):
print '{0} is in {1}:'.format(value, interval), value in interval
With high precision:
In [18]:
n = 1e-100
with mp.workdps(100):
P0 = mp.matrix([[1-3*(mp.exp(n)-1), 3*(mp.exp(n)-1)], [mp.exp(n)-1, 1-(mp.exp(n)-1)]])
print P0
In [19]:
print mp.mp.stoch_eig(P0)
gth_solve returns the (normalized) nontrivial solution to $x A = 0$ for
an irreducible
transition rate martix A.
In [20]:
A = P - mp.eye(3)
In [21]:
print A
In [22]:
y = mp.mp.gth_solve(A)
In [23]:
print y
In [24]:
print mp.chop(y * A)
In [25]:
print mp.norm(y*A, p=1)
In [26]:
print mp.eps
Stochastic matrices arising from a certain game-theoretic model, called the Kandori-Mailath-Rob (KMR) model
In [27]:
from test_eigen_markov import KMR_Markov_matrix_sequential
In [28]:
KMR_matrices = \
[KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14),
KMR_Markov_matrix_sequential(N=10, p=1./3, epsilon=1e-10),
KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)]
In [29]:
print 'KMR_matrices[0] ({0} x {1})'.format(KMR_matrices[0].rows, KMR_matrices[0].cols)
print KMR_matrices[0]
In [30]:
print 'KMR_matrices[1] ({0} x {1})'.format(KMR_matrices[1].rows, KMR_matrices[0].cols)
mp.nprint(KMR_matrices[1], 9)
In [31]:
print 'KMR_matrices[2] ({0} x {1})'.format(KMR_matrices[2].rows, KMR_matrices[2].cols)
In [32]:
x0 = mp.mp.stoch_eig(KMR_matrices[0])
print x0
In [33]:
print x0 * KMR_matrices[0]
In [34]:
print mp.norm(x0*KMR_matrices[0] - x0, p=1)
In [35]:
vecs0 = []
for P in KMR_matrices:
vecs0.append(mp.mp.stoch_eig(P))
Errors $\lVert x P - x\rVert_1$:
In [36]:
for i, P in enumerate(KMR_matrices):
print 'KMR_matrices[{0}]:'.format(i), mp.norm(vecs0[i]*P - vecs0[i], p=1)
Eigenvectors can be found by eig, a routine for the general eigenvalue problem.
Here we are looking for the dominant eigenvector of a stochastic matrix.
In [37]:
def stoch_eig_by_eig(P):
E, EL = mp.eig(mp.matrix(P), left=True, right=False)
E, EL = mp.eig_sort(E, EL)
dom_eigvec = EL[EL.rows-1, :]
return dom_eigvec/sum(dom_eigvec)
In [38]:
P = mp.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]])
print P
In [39]:
z = stoch_eig_by_eig(P)
print z
In [40]:
print z * P
In [41]:
mp.norm(z*P - z, p=1)
Out[41]:
It is natural to suspect that the stoch_eig routine, which is specialized for stochastic matrices,
performs better than the general purpose routine eig for those matrices.
It seems it is indeed the case in terms of speed,
and for accuracy, stoch_eig seems at least as accurate as eig.
In [42]:
import time
In [43]:
v0, v1 = [], []
t0, t1 = [], []
for P in KMR_matrices:
start = time.time()
v = mp.mp.stoch_eig(P)
timediff = time.time() - start
v0.append(v)
t0.append(timediff)
for P in KMR_matrices:
start = time.time()
v = stoch_eig_by_eig(P)
timediff = time.time() - start
v1.append(v)
t1.append(timediff)
In [44]:
for i, P in enumerate(KMR_matrices):
print 'KMR_matrices[{0}] ({1} x {2})'.format(i, P.rows, P.cols)
print ' stoch_eig:', mp.norm(v0[i]*P - v0[i], 1), t0[i]
print ' eig :', mp.norm(v1[i]*P - v1[i], 1), t1[i]
In [45]:
for i, P in enumerate(KMR_matrices):
print 'KMR_matrices[{0}] ({1} x {2})'.format(i, P.rows, P.cols)
%timeit mp.mp.stoch_eig(P)
%timeit stoch_eig_by_eig(P)
In [46]:
sizes = [10, 30, 50]
rand_matrices = []
for n in sizes:
P = mp.randmatrix(n, n)
for i in range(n):
P[i, :] /= mp.fsum(P[i, :])
rand_matrices.append(P)
In [47]:
v0, v1 = [], []
t0, t1 = [], []
for P in rand_matrices:
start = time.time()
v = mp.mp.stoch_eig(P)
timediff = time.time() - start
v0.append(v)
t0.append(timediff)
for P in rand_matrices:
start = time.time()
v = stoch_eig_by_eig(P)
timediff = time.time() - start
v1.append(v)
t1.append(timediff)
In [48]:
for i, P in enumerate(rand_matrices):
print 'rand_matrices[{0}] ({1} x {2})'.format(i, P.rows, P.cols)
print ' stoch_eig:', mp.norm(v0[i]*P - v0[i], 1), t0[i]
print ' eig :', mp.norm(v1[i]*P - v1[i], 1), t1[i]
In [49]:
for i, P in enumerate(rand_matrices):
print 'rand_matrices[{0}] ({1} x {2})'.format(i, P.rows, P.cols)
%timeit mp.mp.stoch_eig(P)
%timeit stoch_eig_by_eig(P)
In [50]:
import os.path
import sys
print("mpmath imported from %s" % os.path.dirname(mp.__file__))
print("mpmath backend: %s" % mp.libmp.backend.BACKEND)
print("mpmath mp class: %s" % repr(mp.mp))
print("mpmath version: %s" % mp.__version__)
print("Python version: %s" % sys.version)
In [50]: