(Using sympy.mpmath)
import numpy as np
from mc_compute_stationary_mpmath import mc_compute_stationary_mpmath
A = [[.4, .6], [.2, .8]]
mc_compute_stationary_mpmath(A) # default precision = 17
[array([mpf('0.25'), mpf('0.75')], dtype=object)]
type(_[0][0])
sympy.mpmath.ctx_mp_python.mpf
mc_compute_stationary_mpmath(A, precision=100)
[array([mpf('0.25'), mpf('0.75')], dtype=object)]
B = np.identity(2)
mc_compute_stationary_mpmath(B)
[array([mpf('1.0'), mpf('0.0')], dtype=object), array([mpf('0.0'), mpf('1.0')], dtype=object)]
C = np.identity(3)
mc_compute_stationary_mpmath(C)
[array([mpf('1.0'), mpf('0.0'), mpf('0.0')], dtype=object), array([mpf('0.0'), mpf('1.0'), mpf('0.0')], dtype=object), array([mpf('0.0'), mpf('0.0'), mpf('1.0')], 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)
eigvecs_P = mc_compute_stationary_mpmath(P)
eigvecs_P
[array([mpf('5.0007905383225959e-15'), mpf('7.5003952691612597e-29'), mpf('1.4999999999999773e-14'), mpf('0.99999999999998')], dtype=object)]
[np.dot(vec, P) for vec in eigvecs_P]
[array([mpf('5.0007905383225958e-15'), mpf('7.5003952691612591e-29'), mpf('1.4999999999999772e-14'), mpf('0.99999999999998002')], dtype=object)]
Q = KMR_Markov_matrix_simultaneous(N=5, p=1./3, epsilon=1e-15)
mc_compute_stationary_mpmath(Q)
[array([mpf('5.2710989716153529e-77'), mpf('4.7477838728799111e-61'), mpf('1.710569414459031e-45'), mpf('3.0814879110196131e-30'), mpf('2.775557561562889e-15'), mpf('0.99999999999999722')], dtype=object)]
mc_compute_stationary_mpmath(Q, precision=100)
[array([mpf('8.8722137051968959e-45'), mpf('2.2655312650280333e-59'), mpf('1.710569414459031e-45'), mpf('3.0814879110196131e-30'), mpf('2.775557561562889e-15'), mpf('0.99999999999999722')], dtype=object)]
R = KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)
eigvecs_R = mc_compute_stationary_mpmath(R)
eigvecs_R
[array([mpf('1.7806708906869897e-21'), mpf('2.4159856305803426e-22'), mpf('1.578282070228375e-23'), mpf('6.6092213996170641e-25'), mpf('1.9927300702419417e-26'), mpf('4.6063107155823758e-28'), mpf('8.4873397144475981e-30'), mpf('1.2794987709652257e-31'), mpf('1.6074386201240665e-33'), mpf('1.707766490924184e-35'), mpf('6.1216883692955772e-33'), mpf('1.8827013498726174e-30'), mpf('4.9954342797262025e-28'), mpf('1.1470285634880882e-25'), mpf('2.282586841341555e-23'), mpf('3.9367014390337377e-21'), mpf('5.8755268977578531e-19'), mpf('7.5656049289364331e-17'), mpf('8.3641965603241664e-15'), mpf('7.8843558102845161e-13'), mpf('6.2759472249864753e-11'), mpf('4.1630449925743605e-9'), mpf('2.259398055060813e-7'), mpf('9.7743524555891698e-6'), mpf('0.00032418268977704092'), mpf('0.0077414826318757366'), mpf('0.11850423413409787'), mpf('0.87342009602538692')], dtype=object)]
[np.dot(vec, R) for vec in eigvecs_R]
[array([mpf('1.7806708906869898e-21'), mpf('2.4159856305803426e-22'), mpf('1.5782820702283751e-23'), mpf('6.6092213996170639e-25'), mpf('1.9927300702419399e-26'), mpf('4.606310715582319e-28'), mpf('8.4873397144518295e-30'), mpf('1.2794987708712456e-31'), mpf('1.6074386243394429e-33'), mpf('1.7077650819974548e-35'), mpf('6.1216880743481326e-33'), mpf('1.8827013488692391e-30'), mpf('4.9954342797235464e-28'), mpf('1.1470285634880883e-25'), mpf('2.2825868413415551e-23'), mpf('3.9367014390337374e-21'), mpf('5.8755268977578527e-19'), mpf('7.5656049289364331e-17'), mpf('8.3641965603241669e-15'), mpf('7.884355810284516e-13'), mpf('6.2759472249864745e-11'), mpf('4.16304499257436e-9'), mpf('2.2593980550608133e-7'), mpf('9.7743524555891707e-6'), mpf('0.00032418268977704094'), mpf('0.0077414826318757358'), mpf('0.11850423413409787'), mpf('0.87342009602538684')], dtype=object)]
import sympy.mpmath as mp
mp.__version__
'0.18'
The working precision is controlled by a context object called mp
.
See: Setting the precision
print mp.mp
Mpmath settings: mp.prec = 53 [default: 53] mp.dps = 15 [default: 15] mp.trap_complex = False [default: False]
The routine eig
solves the eigenvalue problem.
See: The eigenvalue problem
E_A, EL_A = mp.eig(mp.matrix(A), left=True, right=False)
E_A, EL_A[0, :]/sum(EL_A[0, :]), EL_A[1, :]/sum(EL_A[1, :])
([mpf('0.20000000000000001'), mpf('1.0')], matrix( [['7.12081624598818e+15', '-7.12081624598818e+15']]), matrix( [['0.25', '0.75']]))
Let us see what happens if we increase the precision.
See: Setting the precision
mp.mp.dps += 15
E_A, EL_A = mp.eig(mp.matrix(A), left=True, right=False)
E_A, EL_A[0, :]/sum(EL_A[0, :]), EL_A[1, :]/sum(EL_A[1, :])
([mpf('0.200000000000000024980018054066121'), mpf('1.00000000000000004163336342344357')], matrix( [['14411518807585573.1730559644772', '-14411518807585572.1730559644772']]), matrix( [['0.250000000000000004336808689942', '0.749999999999999995663191310058']]))
Arithemtic seems to be more accurate if the mpf
values are created from strings.
See: Providing correct input
AA = np.array(A).astype('str')
AA
array([['0.4', '0.6'], ['0.2', '0.8']], dtype='|S32')
E_AA, EL_AA = mp.eig(mp.matrix(AA), left=True, right=False)
E_AA, EL_AA
([mpf('0.200000000000000000000000000000015'), mpf('1.0')], matrix( [['0.790569415042094832999723386108', '-0.790569415042094832999723386108'], ['0.316227766016837933199889354443', '0.94868329805051379959966806333']]))
EL_AA[1, :]/sum(EL_AA[1, :])
matrix( [['0.25', '0.75']])
KMR Markov matrices:
E_P, EL_P = mp.eig(mp.matrix(P), left=True, right=False)
E_P, EL_P
([mpf('-3.99638918679531821285541014806953e-18'), mpf('0.999999999999996670663055853461674'), mpf('0.66666666666666329765802890809357'), mpf('1.00000000000000000399638918679477')], matrix( [['0.333333333333331091275264305103', '-0.999999999999993333333333333361', '0.999999999999993329336944146565', '-0.333333333333331089943134576171'], ['0.999999999999994999999999999925', '4.99999999999999995413265490448e-15', '-4.99999999999999764809418455726e-15', '-0.999999999999994854238180931096'], ['4.99999999999994929582193330112e-15', '-4.9999999999999500344731847569e-15', '-0.999999999999985000000000000254', '0.999999999999984829470157119688'], ['-5.0000000000001000777919432869e-15', '-7.50000000000011167190087028916e-29', '-1.50000000000000731970634277783e-14', '-0.999999999999999999999999999875']]))
E_Q, EL_Q = mp.eig(mp.matrix(Q), left=True, right=False)
E_Q
[mpf('-1.49690297038509714586703443814896e-46'), mpf('1.49690297038508838775163240798528e-46'), mpf('-8.23548869832210537509734489159652e-47'), mpf('0.999999999999999946487043362145377'), mpf('7.00649232162408535461864791644958e-46'), mpf('1.00000000000000000000000000000079')]
E_R, EL_R = mp.eig(mp.matrix(R), left=True, right=False)
E_R
[mpf('3.98920755789498887221462802905316e-18'), mpf('0.347402974926105442776856364857058'), mpf('0.401430314822682992223009748578764'), mpf('0.44346370835574638375517773575688'), mpf('0.481394103526276933318667259327033'), mpf('0.518514068942550355629297463648163'), mpf('0.555555410535722954392186901330341'), mpf('0.59259258931721140166631097759314'), mpf('0.680085877238326372501428175726829'), mpf('0.629629629575695062410538976338045'), mpf('0.735169997185588140710017801414445'), mpf('0.666666666665970246357253600572777'), mpf('0.777031877313813931845128452635375'), mpf('0.703703703703698035299316901128785'), mpf('0.740740740740740683017994145687565'), mpf('0.814767949255607752872064249682868'), mpf('0.777777777777777774151331208252009'), mpf('0.851850411773480265661651320106536'), mpf('0.814814814814814765566251882406014'), mpf('0.888888865223574719753626902501014'), mpf('0.851851851851851824978215819009429'), mpf('0.925925925713011059448614773704718'), mpf('0.962962962961977024988484907626878'), mpf('0.88888888888888886677964618523252'), mpf('0.925925925925925916778464982765757'), mpf('0.962962962962962979355951480225381'), mpf('0.999999999999998167348029338613544'), mpf('1.00000000000000000015687216788222')]