In [1]:

    

from __future__ import division

import numpy as np
from scipy import sparse
from numba import jit


    

In [2]:

    

@jit(nopython=True)
def searchsorted(a, v):
    lo = -1
    hi = len(a)
    while(lo < hi-1):
        m = (lo + hi) // 2
        if v < a[m]:
            hi = m
        else:
            lo = m
    return hi


    

In [3]:

    

@jit
def mc_sample_path_sp(P, init=0, sample_size=1000):
    P = sparse.csr_matrix(P)
    n = P.shape[0]
    data = P.data
    indices = P.indices
    indptr = P.indptr

    # CDFs, one for each row of P
    cdfs1d = np.empty(P.nnz)
    for i in range(n):
        cdfs1d[indptr[i]:indptr[i+1]] = data[indptr[i]:indptr[i+1]].cumsum()

    # Random values, uniformly sampled from [0, 1)
    u = np.random.random(size=sample_size)

    # === set up array to store output === #
    X = np.empty(sample_size, dtype=int)
    if isinstance(init, int):
        X[0] = init
    else:
        cdf0 = np.cumsum(init)
        X[0] = searchsorted(cdf0, u[0])

    # === generate the sample path === #
    for t in range(sample_size-1):
        k = searchsorted(cdfs1d[indptr[X[t]]:indptr[X[t]+1]], u[t+1])
        X[t+1] = indices[indptr[X[t]]+k]

    return X


    

In [4]:

    

P = np.array([[0.4, 0.6], [0.2, 0.8]])
init = (0.25, 0.75)


    

In [5]:

    

mc_sample_path_sp(P, init=init, sample_size=10)


    

    
Out[5]:





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

In [6]:

    

sample_size = 10**4
X = mc_sample_path_sp(P, init=init, sample_size=sample_size)
print X.sum() / sample_size


    

    




0.7514

In [7]:

    

P = sparse.csr_matrix([[0.4, 0.6], [0.2, 0.8]])
init = (0.25, 0.75)
sample_size = 10**5 * 2
%timeit mc_sample_path_sp(P, init=init, sample_size=sample_size)


    

    




1 loops, best of 3: 4.71 ms per loop

In [8]:

    

import platform
print platform.platform()


    

    




Darwin-13.4.0-x86_64-i386-64bit

In [9]:

    

import sys
print sys.version


    

    




2.7.8 (default, Jul  2 2014, 10:14:46) 
[GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)]

In [10]:

    

print np.__version__


    

    




1.9.0

In [11]:

    

import scipy
print scipy.__version__


    

    




0.14.0

In [12]:

    

import numba
print numba.__version__


    

    




0.18.2

In [12]: