In [ ]:

    

#about DiscreteRV
from numpy import cumsum
from numpy.random import uniform


class DiscreteRV(object):
    """
    Generates an array of draws from a discrete random variable with
    vector of probabilities given by q.
    Parameters
    ----------
    q : array_like(float)
        Nonnegative numbers that sum to 1
    Attributes
    ----------
    q : see Parameters
    Q : array_like(float)
        The cumulative sum of q
    """

    def __init__(self, q):
        self._q = q
        self.Q = cumsum(q)

    def __repr__(self):
        return "DiscreteRV with {n} elements".format(n=self._q.size)

    def __str__(self):
        return self.__repr__()

    @property
    def q(self):
        """
        Getter method for q.
        """
        return self._q

    @q.setter
    def q(self, val):
        """
        Setter method for q.
        """
        self._q = val
        self.Q = cumsum(val)

    def draw(self, k=1):
        """
        Returns k draws from q.
        For each such draw, the value i is returned with probability
        q[i].
        Parameters
        -----------
        k : scalar(int), optional
            Number of draws to be returned
        Returns
        -------
        array_like(int)
            An array of k independent draws from q
        """
        return self.Q.searchsorted(uniform(0, 1, size=k))
    
    #qは各状態（整数で表記）の持つ確率の列。Qはその累積確率を表示した列。(0,1)の一様分布からとってきたいサンプル数だけ数値を引き出す。
    #一様分布からとってきた値を累積確率の列にorderを正しく挿入することを考えたときに、それぞれ何番目に入るのかをdrawで出力する。
    #状態を整数（index）で表記しているので、この時出力されるorderが、規定の確率分布に対して得られた状態となっている。


    

In [2]:

    

from __future__ import division
from quantecon import DiscreteRV 

#class DiscreteRVのobjectは合計が1になるような確率の列（tuple）
psi = (0.1, 0.9)
d = DiscreteRV(psi)

#drawの引数はとってきたいサンプル数
d.draw(5)


    

    
Out[2]:





array([1, 1, 1, 1, 1])

In [3]:

    

# asarrayを試す
import numpy as np
import quantecon as qe

P = [[0.4,0.6],[0.2,0.8]]


    

In [4]:

    

P


    

    
Out[4]:





[[0.4, 0.6], [0.2, 0.8]]

In [5]:

    

P = np.asarray(P)


    

In [6]:

    

P


    

    
Out[6]:





array([[ 0.4,  0.6],
       [ 0.2,  0.8]])

In [8]:

    

P[1,0]


    

    
Out[8]:





0.20000000000000001

In [9]:

    

P[0,1]


    

    
Out[9]:





0.59999999999999998

In [10]:

    

P[1,:]


    

    
Out[10]:





array([ 0.2,  0.8])

In [11]:

    

P[:,0]


    

    
Out[11]:





array([ 0.4,  0.2])

In [12]:

    

len(P)


    

    
Out[12]:





2

In [23]:

    

# make path function

def mc_sample_path(P, init = 0, sample_size = 1000):
    # setting
    P = np.asarray(P)
    n = len(P)
    X = np.empty(sample_size, dtype = int)
    X[0] = init
    
    # Markov Chain
    conditional_dist = [qe.DiscreteRV(P[i,:]) for i in range(n)]
    for t in range(sample_size-1):    # Xのサイズが1000なので、-1しないと1001が入らない
        X[t+1] = conditional_dist[X[t]].draw()
        
    return X


    

In [24]:

    

mc_sample_path(P)


    

    
Out[24]:





array([0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0,
       1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1,
       1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
       1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
       0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0,
       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
       1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
       1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1,
       0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1,
       1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
       1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
       1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
       0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
       1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1,
       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1,
       1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0,
       0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1,
       0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1,
       0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0,
       0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
       1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
       0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1,
       1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0,
       0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1,
       1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1,
       1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0,
       1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0,
       1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
       1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1])

In [25]:

    

X = mc_sample_path(P)
np.mean(X==0)  # X ==0 の割合


    

    
Out[25]:





0.24399999999999999

In [26]:

    

np.mean(X)


    

    
Out[26]:





0.75600000000000001

In [27]:

    

np.mean(X ==1)


    

    
Out[27]:





0.75600000000000001

In [1]:

    

#  quanteconにはすでにMarkov Chainがはいっている
# なにやらできない　

import quantecon as qe
import numpy as np

P = np.asarray([[0.4,0.6],[0.2,0.8]])
q = qe.MarkovChain(P)
Y = q.simulate(ts_length = 10000)
np.mean(Y == 0)


    

    




---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
<ipython-input-1-3e4a0eebda4f> in <module>()
      6 P = np.asarray([[0.4,0.6],[0.2,0.8]])
      7 q = qe.MarkovChain(P)
----> 8 Y = q.simulate(ts_length = 10000)
      9 np.mean(Y == 0)

/Users/susu/anaconda/lib/python2.7/site-packages/quantecon/markov/core.pyc in simulate(self, ts_length, init, num_reps, random_state)
    391         if not self.is_sparse:  # Dense
    392             _generate_sample_paths(
--> 393                 self.cdfs, init_states, random_values, out=X
    394             )
    395         else:  # Sparse

/Users/susu/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in _compile_for_args(self, *args, **kws)
    155         assert not kws
    156         sig = tuple([self.typeof_pyval(a) for a in args])
--> 157         return self.compile(sig)
    158 
    159     def inspect_types(self, file=None):

/Users/susu/anaconda/lib/python2.7/site-packages/numba/dispatcher.pyc in compile(self, sig)
    275                                           self.py_func,
    276                                           args=args, return_type=return_type,
--> 277                                           flags=flags, locals=self.locals)
    278 
    279             # Check typing error if object mode is used

/Users/susu/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(typingctx, targetctx, func, args, return_type, flags, locals, library)
    545     pipeline = Pipeline(typingctx, targetctx, library,
    546                         args, return_type, flags, locals)
--> 547     return pipeline.compile_extra(func)
    548 
    549 

/Users/susu/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_extra(self, func)
    291                 raise e
    292 
--> 293         return self.compile_bytecode(bc, func_attr=self.func_attr)
    294 
    295     def compile_bytecode(self, bc, lifted=(),

/Users/susu/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in compile_bytecode(self, bc, lifted, func_attr)
    299         self.lifted = lifted
    300         self.func_attr = func_attr
--> 301         return self._compile_bytecode()
    302 
    303     def compile_internal(self, bc, func_attr=DEFAULT_FUNCTION_ATTRIBUTES):

/Users/susu/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in _compile_bytecode(self)
    532 
    533         pm.finalize()
--> 534         return pm.run(self.status)
    535 
    536 

/Users/susu/anaconda/lib/python2.7/site-packages/numba/compiler.pyc in run(self, status)
    189                     # No more fallback pipelines?
    190                     if is_final_pipeline:
--> 191                         raise patched_exception
    192                     # Go to next fallback pipeline
    193                     else:

TypingError: Failed at nopython (nopython frontend)
Undeclared getitem(array(float64, 2d, C, nonconst), int64)
File "../../../../../anaconda/lib/python2.7/site-packages/quantecon/markov/core.py", line 438

In [2]:

    

# chack irreducibility

P = [[0.9,0.1,0.0],[0.4,0.4,0.2],[0.1,0.1,0.8]]
mc = qe.MarkovChain(P)
mc.is_irreducible


    

    
Out[2]:





True

In [5]:

    

mc.stationary_distributions


    

    
Out[5]:





array([[ 0.71428571,  0.14285714,  0.14285714]])

In [7]:

    

% matplotlib inline


    

In [15]:

    

# check the convergence to the stationary distribution
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from quantecon import mc_compute_stationary

P = ((0.971, 0.029, 0.000),
     (0.145, 0.778, 0.077),
     (0.000, 0.508, 0.492))
P = np.array(P)

psi = (0.0, 0.8, 0.2)        # Initial condition sumが1になるように変化させても収束する

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_zlim(0, 1)
ax.set_xticks((0.25, 0.5, 0.75))
ax.set_yticks((0.25, 0.5, 0.75))
ax.set_zticks((0.25, 0.5, 0.75))

x_vals, y_vals, z_vals = [], [], []
for t in range(20):
    x_vals.append(psi[0])
    y_vals.append(psi[1])
    z_vals.append(psi[2])
    psi = np.dot(psi, P)

ax.scatter(x_vals, y_vals, z_vals, c='r', s=60)

psi_star = mc_compute_stationary(P)[0]
ax.scatter(psi_star[0], psi_star[1], psi_star[2], c='k', s=60)

plt.show()


    

In [ ]: