mc_compute_stationary_mpmath: demonstration

(Using sympy.mpmath)


In [1]:
import numpy as np
from mc_compute_stationary_mpmath import mc_compute_stationary_mpmath

In [2]:
A = [[.4, .6], [.2, .8]]
mc_compute_stationary_mpmath(A)  # default precision = 17


Out[2]:
[array([mpf('0.25'), mpf('0.75')], dtype=object)]

In [3]:
type(_[0][0])


Out[3]:
sympy.mpmath.ctx_mp_python.mpf

In [4]:
mc_compute_stationary_mpmath(A, precision=100)


Out[4]:
[array([mpf('0.25'), mpf('0.75')], dtype=object)]

In [5]:
B = np.identity(2)
mc_compute_stationary_mpmath(B)


Out[5]:
[array([mpf('1.0'), mpf('0.0')], dtype=object),
 array([mpf('0.0'), mpf('1.0')], dtype=object)]

In [6]:
C = np.identity(3)
mc_compute_stationary_mpmath(C)


Out[6]:
[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)]

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 [7]:
from test_mc_compute_stationary import KMR_Markov_matrix_sequential, KMR_Markov_matrix_simultaneous

In [8]:
P = KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14)
eigvecs_P = mc_compute_stationary_mpmath(P)

In [9]:
eigvecs_P


Out[9]:
[array([mpf('5.0007905383225959e-15'), mpf('7.5003952691612597e-29'),
        mpf('1.4999999999999773e-14'), mpf('0.99999999999998')], dtype=object)]

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


Out[10]:
[array([mpf('5.0007905383225958e-15'), mpf('7.5003952691612591e-29'),
        mpf('1.4999999999999772e-14'), mpf('0.99999999999998002')], dtype=object)]

In [11]:
Q = KMR_Markov_matrix_simultaneous(N=5, p=1./3, epsilon=1e-15)
mc_compute_stationary_mpmath(Q)


Out[11]:
[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)]

In [12]:
mc_compute_stationary_mpmath(Q, precision=100)


Out[12]:
[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)]

In [13]:
R = KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)
eigvecs_R = mc_compute_stationary_mpmath(R)

In [14]:
eigvecs_R


Out[14]:
[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)]

In [15]:
[np.dot(vec, R) for vec in eigvecs_R]


Out[15]:
[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)]

mpmath

mc_compute_stationary_mpmath uses mpmath.eig from mpmath.


In [16]:
import sympy.mpmath as mp

In [17]:
mp.__version__


Out[17]:
'0.18'

The working precision is controlled by a context object called mp.
See: Setting the precision


In [18]:
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


In [19]:
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, :])


Out[19]:
([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


In [20]:
mp.mp.dps += 15

In [21]:
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, :])


Out[21]:
([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


In [22]:
AA = np.array(A).astype('str')
AA


Out[22]:
array([['0.4', '0.6'],
       ['0.2', '0.8']], 
      dtype='|S32')

In [23]:
E_AA, EL_AA = mp.eig(mp.matrix(AA), left=True, right=False)
E_AA, EL_AA


Out[23]:
([mpf('0.200000000000000000000000000000015'), mpf('1.0')], matrix(
 [['0.790569415042094832999723386108', '-0.790569415042094832999723386108'],
  ['0.316227766016837933199889354443', '0.94868329805051379959966806333']]))

In [24]:
EL_AA[1, :]/sum(EL_AA[1, :])


Out[24]:
matrix(
[['0.25', '0.75']])

KMR Markov matrices:


In [25]:
E_P, EL_P = mp.eig(mp.matrix(P), left=True, right=False)

In [26]:
E_P, EL_P


Out[26]:
([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']]))

In [27]:
E_Q, EL_Q = mp.eig(mp.matrix(Q), left=True, right=False)

In [28]:
E_Q


Out[28]:
[mpf('-1.49690297038509714586703443814896e-46'),
 mpf('1.49690297038508838775163240798528e-46'),
 mpf('-8.23548869832210537509734489159652e-47'),
 mpf('0.999999999999999946487043362145377'),
 mpf('7.00649232162408535461864791644958e-46'),
 mpf('1.00000000000000000000000000000079')]

In [29]:
E_R, EL_R = mp.eig(mp.matrix(R), left=True, right=False)

In [30]:
E_R


Out[30]:
[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')]

In [30]: