Generates one sample path from the Markov chain represented by
    (n x n) transition matrix P on state space S = {{0,...,n-1}}.

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

    init : array_like(float ndim=1) or scalar(int), optional(default=0)
        If init is an array_like, then it is treated as the initial
        distribution across states.  If init is a scalar, then it
        treated as the deterministic initial state.

    sample_size : scalar(int), optional(default=1000)
        The length of the sample path.

    random_state : int or np.random.RandomState, optional
        Random seed (integer) or np.random.RandomState instance to set
        the initial state of the random number generator for
        reproducibility. If None, a randomly initialized RandomState is
        used.

    Returns
    -------
    X : array_like(int, ndim=1)
        The simulation of states.

def mc_sample_path(P, init=0, sample_size=1000, random_state=None):
    """
    Generates one sample path from the Markov chain represented by
    (n x n) transition matrix P on state space S = {{0,...,n-1}}.

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

    init : array_like(float ndim=1) or scalar(int), optional(default=0)
        If init is an array_like, then it is treated as the initial
        distribution across states.  If init is a scalar, then it
        treated as the deterministic initial state.

    sample_size : scalar(int), optional(default=1000)
        The length of the sample path.

    random_state : int or np.random.RandomState, optional
        Random seed (integer) or np.random.RandomState instance to set
        the initial state of the random number generator for
        reproducibility. If None, a randomly initialized RandomState is
        used.

    Returns
    -------
    X : array_like(int, ndim=1)
        The simulation of states.

    """
    random_state = check_random_state(random_state)

    if isinstance(init, numbers.Integral):
        X_0 = init
    else:
        cdf0 = np.cumsum(init)
        u_0 = random_state.random_sample()
        X_0 = searchsorted(cdf0, u_0)

    mc = MarkovChain(P)
    return mc.simulate(ts_length=sample_size, init=X_0,
                       random_state=random_state)