mc_compute_stationary_sympy: demonstration



In [1]:

    
import numpy as np
from mc_compute_stationary_sympy import mc_compute_stationary_sympy



In [2]:

    
A = [[.4, .6], [.2, .8]]
mc_compute_stationary_sympy(A)









    Out[2]:





([1.00000000000000],
 [array([0.250000000000000, 0.750000000000000], dtype=object)])



In [3]:

    
type(_[1][0][0])









    Out[3]:





sympy.core.numbers.Float



In [4]:

    
B = np.identity(2)
mc_compute_stationary_sympy(B)









    Out[4]:





([1.00000000000000],
 [array([1.00000000000000, 0], dtype=object),
  array([0, 1.00000000000000], dtype=object)])



In [5]:

    
C = np.identity(3)
mc_compute_stationary_sympy(C)









    Out[5]:





([1.00000000000000],
 [array([1.00000000000000, 0, 0], dtype=object),
  array([0, 1.00000000000000, 0], dtype=object),
  array([0, 0, 1.00000000000000], dtype=object)])

KMR Markov matrices

Let us try matrices generated by the KMR model where action 1 is risk-dominant, so that the theory says that the unique stationary distribution is close to the distribution that concentrates on the state $N$ (where the state is the number of players playing action 1).



In [6]:

    
from test_mc_compute_stationary import KMR_Markov_matrix_sequential, KMR_Markov_matrix_simultaneous



In [7]:

    
P = KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14)
eigvals_P, eigvecs_P = mc_compute_stationary_sympy(P)



In [8]:

    
eigvals_P, eigvecs_P









    Out[8]:





([1.00000000000000, 1.00000000000000],
 [array([-428914250225761., -2.14457125112880, 2.14457125112881,
         428914250225762.], dtype=object),
  array([4.99999999999991e-15, 7.49999999999982e-29, 1.49999999999996e-14,
         0.999999999999980], dtype=object)])



In [9]:

    
[val.evalf(50) for val in eigvals_P]









    Out[9]:





[1.0000000000000017763568394002504646778106689453125,
 1.0000000000000048849813083506887778639793395996094]



In [10]:

    
[np.dot(vec, P) for vec in eigvecs_P]









    Out[10]:





[array([-428914250225759., -2.14457125112881, 2.14457125112881,
        428914250225760.], dtype=object),
 array([4.99999999999991e-15, 7.49999999999987e-29, 1.49999999999997e-14,
        0.999999999999980], dtype=object)]



In [11]:

    
Q = KMR_Markov_matrix_simultaneous(N=5, p=1./3, epsilon=1e-15)
eigvals_Q, eigvecs_Q = mc_compute_stationary_sympy(Q)



In [12]:

    
eigvals_Q, eigvecs_Q









    Out[12]:





([1.00000000000000, 1.00000000000000],
 [array([1.00000000000001, 2.50000000000002e-15, 2.50000000000003e-30,
         -2.79568440577553e-44, -2.51812864961921e-29, -9.07251459847683e-15], dtype=object),
  array([0, 0, 0, 3.08148791101960e-30, 2.77555756156288e-15,
         0.999999999999997], dtype=object)])



In [13]:

    
[val.evalf(50) for val in eigvals_Q]









    Out[13]:





[1.0000000000000024424906541753443889319896697998047,
 1.0000000000000026645352591003756970167160034179688]



In [14]:

    
[np.dot(vec, Q) for vec in eigvecs_Q]









    Out[14]:





[array([1.00000000000001, 2.50000000000003e-15, 2.50000000000003e-30,
        -2.67068440577554e-44, -2.51812864961921e-29, -9.07251459847683e-15], dtype=object),
 array([5.27109897161535e-77, 4.74778387287991e-61, 1.71056941445903e-45,
        3.08148791101961e-30, 2.77555756156289e-15, 0.999999999999997], dtype=object)]



In [15]:

    
R = KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)
eigvals_R, eigvecs_R = mc_compute_stationary_sympy(R)



In [16]:

    
eigvals_R, eigvecs_R









    Out[16]:





([], [])

mc_compute_stationary_sympy fails to return any eigenvalue or eigenvector for this $28 \times 28$ matrix.

SymPy

mc_compute_stationary_sympy uses Matrix.eigenvects() from SymPy.



In [17]:

    
from sympy.matrices import Matrix



In [18]:

    
Matrix(P).transpose().eigenvects()









    Out[18]:





[(1.00000000000000, 1, [Matrix([
   [   -0.999999999999998],
   [-4.99999999999999e-15],
   [              5.0e-15],
   [                  1.0]])]), (1.00000000000000, 1, [Matrix([
   [5.00000000000001e-15],
   [7.49999999999997e-29],
   [1.49999999999999e-14],
   [                 1.0]])]), (0.666666666666665, 1, [Matrix([
   [ 5.00000000000001e-15],
   [-5.00000000000003e-15],
   [                 -1.0],
   [                  1.0]])]), (-3.99638918680497e-18, 1, [Matrix([
   [-0.999999999999998],
   [  2.99999999999999],
   [              -3.0],
   [               1.0]])])]



In [19]:

    
Matrix(Q).transpose().eigenvects()









    Out[19]:





[(0, 4, [Matrix([
   [-1.0],
   [ 1.0],
   [   0],
   [   0],
   [   0],
   [   0]]), Matrix([
   [   0],
   [   0],
   [-1.0],
   [ 1.0],
   [   0],
   [   0]]), Matrix([
   [   0],
   [   0],
   [-1.0],
   [   0],
   [ 1.0],
   [   0]]), Matrix([
   [   0],
   [   0],
   [-1.0],
   [   0],
   [   0],
   [ 1.0]])]), (1.00000000000000, 1, [Matrix([
   [   -110223024625157.0],
   [   -0.275557561562892],
   [-2.75557561562893e-16],
   [ 3.08148791101961e-30],
   [ 2.77555756156289e-15],
   [                  1.0]])]), (1.00000000000000, 1, [Matrix([
   [                   0],
   [                   0],
   [                   0],
   [3.08148791101961e-30],
   [2.77555756156289e-15],
   [                 1.0]])])]



In [20]:

    
Matrix(R).transpose().eigenvects()









    Out[20]:





[]



In [20]: