import numpy as np
from mc_compute_stationary_sympy import mc_compute_stationary_sympy
A = [[.4, .6], [.2, .8]]
mc_compute_stationary_sympy(A)
([1.00000000000000], [array([0.250000000000000, 0.750000000000000], dtype=object)])
type(_[1][0][0])
sympy.core.numbers.Float
B = np.identity(2)
mc_compute_stationary_sympy(B)
([1.00000000000000], [array([1.00000000000000, 0], dtype=object), array([0, 1.00000000000000], dtype=object)])
C = np.identity(3)
mc_compute_stationary_sympy(C)
([1.00000000000000], [array([1.00000000000000, 0, 0], dtype=object), array([0, 1.00000000000000, 0], dtype=object), array([0, 0, 1.00000000000000], dtype=object)])
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).
from test_mc_compute_stationary import KMR_Markov_matrix_sequential, KMR_Markov_matrix_simultaneous
P = KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14)
eigvals_P, eigvecs_P = mc_compute_stationary_sympy(P)
eigvals_P, eigvecs_P
([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)])
[val.evalf(50) for val in eigvals_P]
[1.0000000000000017763568394002504646778106689453125, 1.0000000000000048849813083506887778639793395996094]
[np.dot(vec, P) for vec in eigvecs_P]
[array([-428914250225759., -2.14457125112881, 2.14457125112881, 428914250225760.], dtype=object), array([4.99999999999991e-15, 7.49999999999987e-29, 1.49999999999997e-14, 0.999999999999980], dtype=object)]
Q = KMR_Markov_matrix_simultaneous(N=5, p=1./3, epsilon=1e-15)
eigvals_Q, eigvecs_Q = mc_compute_stationary_sympy(Q)
eigvals_Q, eigvecs_Q
([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)])
[val.evalf(50) for val in eigvals_Q]
[1.0000000000000024424906541753443889319896697998047, 1.0000000000000026645352591003756970167160034179688]
[np.dot(vec, Q) for vec in eigvecs_Q]
[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)]
R = KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)
eigvals_R, eigvecs_R = mc_compute_stationary_sympy(R)
eigvals_R, eigvecs_R
([], [])
mc_compute_stationary_sympy fails to return any eigenvalue or eigenvector for this $28 \times 28$ matrix.
mc_compute_stationary_sympy uses Matrix.eigenvects()
from SymPy.
from sympy.matrices import Matrix
Matrix(P).transpose().eigenvects()
[(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]])])]
Matrix(Q).transpose().eigenvects()
[(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]])])]
Matrix(R).transpose().eigenvects()
[]