In [1]:

    

"""
Filename: mc_compute_stationary_mpmath.py

Author: Daisuke Oyama

"""
import sys
import numpy as np

try:
    try:
        import mpmath as mp
    except ImportError:
        import sympy.mpmath as mp
except ImportError:
    sys.stderr.write('Failed to import mpmath\n')
    sys.exit(1)


def mc_compute_stationary_mpmath(P, precision=17, irreducible=False, ltol=0, utol=None):
    """
    Computes the stationary distributions of Markov matrix P.

    Parameters
    ----------
    P : array_like(float, ndim=2)
        A discrete Markov transition matrix.

    precision : scalar(int), optional(default: 17)
        Decimal precision in float-point arithemetic with mpmath.
        mpmath.mp.dps is set to *precision*.

    irreducible : bool, optional(default: False)
        Set True if P is known a priori to be irreducible
        (for any i, j, (P^k)_{ij} > 0 for some k).
        If True, the eigenvector for the maximum eigenvalue is returned.

    ltol, utol: scalar(float), optional(default: ltol=0, utol=None)
        Lower and upper tolerance levels.
        Find eigenvectors for eigenvalues in [1-ltol, 1+utol]
        (where [1-ltol, 1+utol] = [1-ltol, +inf) when utol=None).

    Returns
    -------
    vecs : list of numpy.arrays of mpmath.ctx_mp_python.mpf
        A list of the eigenvectors of whose eigenvalues in [1-ltol, 1+utol].

    Notes
    -----
    mpmath 0.18 or above is required.

    References
    ----------

        http://mpmath.org/doc/current

    """
    LTOL = ltol  # Lower tolerance level
    if utol is None:  # Upper tolerance level
        UTOL = 'inf'
    else:
        UTOL = utol

    with mp.workdps(precision):  # Temporarily change the working precision
        E, EL = mp.eig(mp.matrix(P), left=True, right=False)
        # E  : a list of length n containing the eigenvalues of A
        # EL : a matrix whose rows contain the left eigenvectors of A
        # See: github.com/fredrik-johansson/mpmath/blob/master/mpmath/matrices/eigen.py
        E, EL = mp.eig_sort(E, EL)  # Sorted in a descending order

        if irreducible:
            num_eigval_one = 1
        else:
            num_eigval_one = sum(
                mp.mpf(1) - mp.mpf(LTOL) <= val <= mp.mpf(1) + mp.mpf(UTOL)
                for val in E
                )

        vecs = [np.array((EL[i, :]/sum(EL[i, :])).tolist()[0])
                for i in range(EL.rows-num_eigval_one, EL.rows)]

    return vecs


    

In [2]:

    

from __future__ import division
import numpy as np


def set_tran_matrix(N=6,p=1/3,epsilon=0.01):  #setting up the transition matrix A
	"""
	N :the number of players
	p : the turning point as to which equilibrium is prefered
	epsilon : tthe possibility that a man becomes "irrational"

	if the ratio of players playing 0 exceeds p, we end up the equilibrium(0,0) in an unperturbed game
	I think it's better if I can compute the value of p in this program.For now,we assume it's exogenous.
	"""
	A = np.zeros((N+1,N+1))
	for n in range(1, N):
		A[n, n-1] = \
	    	(n/N) * (epsilon * (1/2) +
                     (1 - epsilon) * (((n)/(N) < p) + ((n)/(N) == p) * (1/2))
                     )
		A[n, n+1] = \
            ((N-n)/N) * (epsilon * (1/2) +
                         (1 - epsilon) * ((n/(N) > p) + (n/(N) == p) * (1/2))
                         )
		A[n,n] = 1- A[n,n-1] - A[n,n+1]

	A[0,0] = 1 - epsilon/2
	A[0,1] = epsilon/2

	A[N,N-1] = epsilon/2
	A[N,N] = 1 - epsilon/2
	
	return A


    

In [3]:

    

P = set_tran_matrix()


    

In [4]:

    

mc_compute_stationary_mpmath(P)


    

    
Out[4]:





[]

In [5]:

    

P = set_tran_matrix(N=26)


    

In [6]:

    

mc_compute_stationary_mpmath(P)


    

    
Out[6]:





[array([mpf('1.7946972729275572e-21'), mpf('2.3448306078450533e-22'),
        mpf('1.4728835476413624e-23'), mpf('5.9211398900153081e-25'),
        mpf('1.7108821290228975e-26'), mpf('3.7828549585157268e-28'),
        mpf('6.6532624792810135e-30'), mpf('9.5524218977083442e-32'),
        mpf('1.1400424773350582e-33'), mpf('2.2786936807187899e-33'),
        mpf('7.7088088919014376e-31'), mpf('2.2313497621135247e-28'),
        mpf('5.55048253316502e-26'), mpf('1.1895111028766731e-23'),
        mpf('2.1980465879585375e-21'), mpf('3.4992901680299916e-19'),
        mpf('4.7874663611360307e-17'), mpf('5.604151799212178e-15'),
        mpf('5.5761310402161177e-13'), mpf('4.6722108505389789e-11'),
        mpf('3.2541948574003983e-9'), mpf('1.850242218921941e-7'),
        mpf('8.3681409446696875e-6'), mpf('0.00028961044312856845'),
        mpf('0.0072040597728231406'), mpf('0.11468863158334436'),
        mpf('0.87780914173405714')], dtype=object)]

In [7]:

    

P


    

    
Out[7]:





array([[ 0.995     ,  0.005     ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.03826923,  0.95692308,  0.00480769,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.07653846,  0.91884615,  0.00461538,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.11480769,  0.88076923,  0.00442308,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.15307692,  0.84269231,
         0.00423077,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.19134615,
         0.80461538,  0.00403846,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.22961538,  0.76653846,  0.00384615,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.26788462,  0.72846154,  0.00365385,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.30615385,  0.69038462,  0.00346154,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.00173077,  0.34769231,
         0.65057692,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.00192308,
         0.38576923,  0.61230769,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.00211538,  0.42384615,  0.57403846,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.00230769,  0.46192308,  0.53576923,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.0025    ,  0.5       ,  0.4975    ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.00269231,  0.53807692,
         0.45923077,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.00288462,
         0.57615385,  0.42096154,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.00307692,  0.61423077,  0.38269231,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.00326923,  0.65230769,  0.34442308,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.00346154,  0.69038462,  0.30615385,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.00365385,  0.72846154,
         0.26788462,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.00384615,
         0.76653846,  0.22961538,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.00403846,  0.80461538,  0.19134615,  0.        ,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.00423077,  0.84269231,  0.15307692,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.00442308,  0.88076923,  0.11480769,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.00461538,  0.91884615,
         0.07653846,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.00480769,
         0.95692308,  0.03826923],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.005     ,  0.995     ]])

In [8]:

    

P = set_tran_matrix()


    

In [9]:

    

mc_compute_stationary_mpmath(P,precision=11)


    

    
Out[9]:





[]

In [10]:

    

mc_compute_stationary_mpmath(P,precision=8)


    

    
Out[10]:





[]

In [11]:

    

mc_compute_stationary_mpmath(P,precision=13)


    

    
Out[11]:





[]

In [ ]: