In [1]:

    

%load mc_sample_path_sp.py


    

In [ ]:

    

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


@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


@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


if __name__ == '__main__':
    P = sparse.csr_matrix([[0.4, 0.6], [0.2, 0.8]])
    init = (0.25, 0.75)
    sample_size = 10
    mc_sample_path_sp(P, init=init, sample_size=sample_size)


    

In [2]:

    

!numba --annotate-html mc_sample_path_sp_annotate.html mc_sample_path_sp.py


    

In [3]:

    

from IPython.display import HTML
HTML('<iframe src=mc_sample_path_sp_annotate.html width=100% height=700></iframe>')


    

    
Out[3]:

In [4]:

    

import platform
print platform.platform()


    

    




Darwin-13.4.0-x86_64-i386-64bit

In [5]:

    

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

    

import numpy
print numpy.__version__


    

    




1.9.0

In [7]:

    

import scipy
print scipy.__version__


    

    




0.14.0

In [8]:

    

import numba
print numba.__version__


    

    




0.18.2

In [8]: