From 4bf577e1ce4068ae5cf124f31cc44c8b056fb682 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 9 May 2015 23:33:05 +0900
Subject: [PATCH 01/12] MarkovChain.simulate API changed

---
 quantecon/markov/core.py            | 201 ++++++++++++++++++----------
 quantecon/markov/tests/test_core.py |  12 +-
 2 files changed, 135 insertions(+), 78 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 22c207389..24e429331 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -140,11 +140,6 @@ class MarkovChain(object):
         List of numpy arrays containing the cyclic classes. Defined only
         when the Markov chain is irreducible.
 
-    Methods
-    -------
-    simulate : Simulates the markov chain for a given initial state or
-        distribution.
-
     """
 
     def __init__(self, P):
@@ -170,6 +165,7 @@ def __init__(self, P):
         self.digraph = DiGraph(self.P)
 
         self._stationary_dists = None
+        self._cdfs = None
 
     def __repr__(self):
         msg = "Markov chain with transition matrix \nP = \n{0}"
@@ -256,98 +252,159 @@ def stationary_distributions(self):
             self._compute_stationary()
         return self._stationary_dists
 
-    def simulate(self, init=0, sample_size=1000):
-        X = mc_sample_path(self.P, init, sample_size)
-        return X
+    @property
+    def cdfs(self):
+        if self._cdfs is None:
+            # See issue #137#issuecomment-96128186
+            cdfs = np.empty((self.n, self.n), order='C')
+            np.cumsum(self.P, axis=-1, out=cdfs)
+            self._cdfs = cdfs
+        return self._cdfs
+
+    def simulate(self, ts_length, init=None, num_reps=None):
+        """
+        Simulate time series of state transitions.
+
+        Parameters
+        ----------
+        ts_length : scalar(int)
+            Length of each simulation.
+
+        init : scalar(int) or array_like(int, ndim=1),
+               optional(default=None)
+            Initial state(s). If None, the initial state is randomly
+            drawn.
+
+        num_reps : scalar(int), optional(default=None)
+            Number of simulations. Relevant only when init is a scalar
+            or None.
+
+        Returns
+        -------
+        X : ndarray(int, ndim=1 or 2)
+            Array containing the sample path(s), of shape (ts_length,)
+            if init is a scalar (integer) or None and num_reps is None;
+            of shape (k, ts_length) otherwise, where k = len(init) if
+            init is an array_like, otherwise k = num_reps.
 
+        """
+        try:
+            num_reps = len(init)  # init is an array; make num_reps not None
+            init_states = np.asarray(init, dtype=int)
+        except:  # init is a scalar(int) or None
+            k = 1 if num_reps is None else num_reps
+            if init is None:
+                init_states = np.random.randint(self.n, size=k)
+            elif isinstance(init, int):
+                init_states = np.ones(k, dtype=int) * init
+            else:
+                raise ValueError(
+                    'init must be int, array_like of ints, or None'
+                )
+
+        X = _simulate_markov_chain(self.cdfs, ts_length, init_states,
+                                   random_state=np.random.RandomState())
+
+        if num_reps is None:
+            return X[0]
+        else:
+            return X
 
-def mc_compute_stationary(P):
+
+def _simulate_markov_chain(P_cdfs, ts_length, init_states, random_state):
     """
-    Computes stationary distributions of P, one for each recurrent
-    class. Any stationary distribution is written as a convex
-    combination of these distributions.
+    Main body of MarkovChain.simulate.
+
+    Generate len(init_states) sample paths of length ts_length.
+
+    Parameters
+    ----------
+    P_cdfs : ndarray(float, ndim=2)
+        Array containing as rows the CDFs of the state transition.
+
+    ts_length : scalar(int)
+        Length of each sample path.
+
+    init_states : array_like(int, ndim=1)
+        Array containing the initial states.
+
+    random_state : np.random.RandomState
 
     Returns
     -------
-    stationary_dists : array_like(float, ndim=2)
-        Array containing the stationary distributions as its rows.
-
-    """
-    return MarkovChain(P).stationary_distributions
+    X : ndarray(int, ndim=2)
+        Array containing the sample paths, of shape (len(init_states),
+        ts_length)
 
+    Notes
+    -----
+    This routine is jit-complied if the module Numba is vailable.
 
-def mc_sample_path(P, init=0, sample_size=1000):
-    """
-    See Section: DocStrings below
     """
-    n = len(P)
+    num_reps = len(init_states)
 
-    # CDFs, one for each row of P
-    cdfs = np.empty((n, n), order='C')  # see issue #137#issuecomment-96128186
-    np.cumsum(P, axis=-1, out=cdfs)
+    # === set up array to store output === #
+    X = np.empty((num_reps, ts_length), dtype=int)
+    X[:, 0] = init_states
 
     # 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])
+    u = random_state.random_sample(size=(num_reps, ts_length-1))
 
-    # === generate the sample path === #
-    for t in range(sample_size-1):
-        X[t+1] = searchsorted(cdfs[X[t]], u[t+1])
+    # === generate the sample paths === #
+    for i in range(num_reps):
+        for t in range(ts_length-1):
+            X[i, t+1] = searchsorted(P_cdfs[X[i, t]], u[i, t])
 
     return X
 
 if numba_installed:
-    mc_sample_path = jit(mc_sample_path)
+    _simulate_markov_chain = jit(_simulate_markov_chain)
+
 
+def mc_compute_stationary(P):
+    """
+    Computes stationary distributions of P, one for each recurrent
+    class. Any stationary distribution is written as a convex
+    combination of these distributions.
 
-# ------------ #
-# -DocStrings- #
-# ------------ #
+    Returns
+    -------
+    stationary_dists : array_like(float, ndim=2)
+        Array containing the stationary distributions as its rows.
 
-# -mc_sample_path() function and MarkovChain.simulate() method- #
-_sample_path_docstr = \
-"""
-Generates one sample path from the Markov chain represented by (n x n)
-transition matrix P on state space S = {{0,...,n-1}}.
+    """
+    return MarkovChain(P).stationary_distributions
 
-Parameters
-----------
-{p_arg}
-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.
+def mc_sample_path(P, init=0, sample_size=1000):
+    """
+    Generates one sample path from the Markov chain represented by
+    (n x n) transition matrix P on state space S = {{0,...,n-1}}.
 
-Returns
--------
-X : array_like(int, ndim=1)
-    The simulation of states.
+    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.
 
-# -Functions- #
+    sample_size : scalar(int), optional(default=1000)
+        The length of the sample path.
 
-# -mc_sample_path- #
-mc_sample_path.__doc__ = _sample_path_docstr.format(p_arg="""
-P : array_like(float, ndim=2)
-    A Markov transition matrix.
-""")
+    Returns
+    -------
+    X : array_like(int, ndim=1)
+        The simulation of states.
 
-# -Methods- #
+    """
+    if isinstance(init, int):
+        X_0 = init
+    else:
+        cdf0 = np.cumsum(init)
+        u_0 = np.random.random(size=1)
+        X_0 = searchsorted(cdf0, u_0)
 
-# -Markovchain.simulate()- #
-if sys.version_info[0] == 3:
-    MarkovChain.simulate.__doc__ = _sample_path_docstr.format(p_arg="")
-elif sys.version_info[0] == 2:
-    MarkovChain.simulate.__func__.__doc__ = \
-        _sample_path_docstr.format(p_arg="")
+    return MarkovChain(P).simulate(ts_length=sample_size, init=X_0)
diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index ace483989..84cab6634 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -226,15 +226,15 @@ def test_simulate_for_matrices_with_C_F_orders():
     P_C = np.array([[0.5, 0.5], [0, 1]], order='C')
     P_F = np.array([[0.5, 0.5], [0, 1]], order='F')
     init = 1
-    sample_size = 10
-    sample_path = np.ones(sample_size, dtype=int)
+    ts_length = 10
+    sample_path = np.ones(ts_length, dtype=int)
 
     computed_C_and_F = \
-        MarkovChain(np.array([[1.]])).simulate(init=0, sample_size=sample_size)
-    assert_array_equal(computed_C_and_F, np.zeros(sample_size, dtype=int))
+        MarkovChain(np.array([[1.]])).simulate(ts_length, init=0)
+    assert_array_equal(computed_C_and_F, np.zeros(ts_length, dtype=int))
 
-    computed_C = MarkovChain(P_C).simulate(init, sample_size)
-    computed_F = MarkovChain(P_F).simulate(init, sample_size)
+    computed_C = MarkovChain(P_C).simulate(ts_length, init)
+    computed_F = MarkovChain(P_F).simulate(ts_length, init=init)
     assert_array_equal(computed_C, sample_path)
     assert_array_equal(computed_F, sample_path)
 

From 4d35930d9c3cdfebfd394db410d99024f9a41704 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 6 Jun 2015 14:41:48 +0900
Subject: [PATCH 02/12] Creating a DiGraph instance delayed

---
 quantecon/markov/core.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 24e429331..27b401f74 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -162,7 +162,7 @@ def __init__(self, P):
         self.n = self.P.shape[0]
 
         # To analyze the structure of P as a directed graph
-        self.digraph = DiGraph(self.P)
+        self._digraph = None
 
         self._stationary_dists = None
         self._cdfs = None
@@ -179,6 +179,12 @@ def __repr__(self):
     def __str__(self):
         return str(self.__repr__)
 
+    @property
+    def digraph(self):
+        if self._digraph is None:
+            self._digraph = DiGraph(self.P)
+        return self._digraph
+
     @property
     def is_irreducible(self):
         return self.digraph.is_strongly_connected

From 4f2ff57b6b22b02530002a74c12b89ed1322a0d1 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 6 Jun 2015 15:29:48 +0900
Subject: [PATCH 03/12] Test for MarkovChain.simulate added

Only for shape
---
 quantecon/markov/tests/test_core.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index 84cab6634..f8688c7fa 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -216,6 +216,23 @@ def test_left_eigen_vec(self):
                 assert_allclose(np.dot(curr_v, mc.P), curr_v, atol=self.TOL)
 
 
+def test_simulate_shape():
+    P = [[0.4, 0.6], [0.2, 0.8]]
+    mc = MarkovChain(P)
+
+    (ts_length, init, num_reps) = (10, None, None)
+    assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
+                       (ts_length,))
+
+    for (ts_length, init, num_reps) in [(10, [0, 1], None), (10, [0, 1], 3)]:
+        assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
+                           (len(init), ts_length))
+
+    for (ts_length, init, num_reps) in [(10, None, 3), (10, None, 1)]:
+        assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
+                           (num_reps, ts_length))
+
+
 def test_simulate_for_matrices_with_C_F_orders():
     """
     Test MarkovChasin.simulate for matrices with C- and F-orders

From 958c7133a2cdc78a58ddc4b79e85afef3c2ba927 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 6 Jun 2015 15:38:21 +0900
Subject: [PATCH 04/12] `random_state` option, as well as tests, added

This will close #148 and close #152
---
 quantecon/markov/core.py            | 42 ++++++++++++++++++++++++-----
 quantecon/markov/tests/test_core.py | 21 +++++++++++++++
 2 files changed, 56 insertions(+), 7 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 27b401f74..a82528375 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -2,7 +2,7 @@
 Authors: Chase Coleman, Spencer Lyon, Daisuke Oyama, Tom Sargent,
          John Stachurski
 
-Filename: core.py
+Filename: mc_tools.py
 
 This file contains some useful objects for handling a finite-state
 discrete-time Markov chain.
@@ -88,11 +88,17 @@
 import numpy as np
 from fractions import gcd
 import sys
+from .graph_tools import DiGraph
 from .gth_solve import gth_solve
-from ..graph_tools import DiGraph
 
 # -Check if Numba is Available- #
+<<<<<<< HEAD
 from ..util import searchsorted, numba_installed, jit
+=======
+from .external import numba_installed, jit
+
+from .util import searchsorted, check_random_state
+>>>>>>> `random_state` option, as well as tests, added
 
 
 class MarkovChain(object):
@@ -107,6 +113,13 @@ class MarkovChain(object):
     P : array_like(float, ndim=2)
         The transition matrix.  Must be of shape n x n.
 
+    random_state : scalar(int) or np.random.RandomState,
+                   optional(default=None)
+        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.
+
     Attributes
     ----------
     P : see Parameters
@@ -142,7 +155,7 @@ class MarkovChain(object):
 
     """
 
-    def __init__(self, P):
+    def __init__(self, P, random_state=None):
         self.P = np.asarray(P)
 
         # Check Properties
@@ -161,6 +174,9 @@ def __init__(self, P):
         # The number of states
         self.n = self.P.shape[0]
 
+        # random state
+        self.random_state = random_state
+
         # To analyze the structure of P as a directed graph
         self._digraph = None
 
@@ -294,6 +310,8 @@ def simulate(self, ts_length, init=None, num_reps=None):
             init is an array_like, otherwise k = num_reps.
 
         """
+        random_state = check_random_state(self.random_state)
+
         try:
             num_reps = len(init)  # init is an array; make num_reps not None
             init_states = np.asarray(init, dtype=int)
@@ -309,7 +327,7 @@ def simulate(self, ts_length, init=None, num_reps=None):
                 )
 
         X = _simulate_markov_chain(self.cdfs, ts_length, init_states,
-                                   random_state=np.random.RandomState())
+                                   random_state=random_state)
 
         if num_reps is None:
             return X[0]
@@ -382,7 +400,7 @@ def mc_compute_stationary(P):
     return MarkovChain(P).stationary_distributions
 
 
-def mc_sample_path(P, init=0, sample_size=1000):
+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}}.
@@ -400,17 +418,27 @@ def mc_sample_path(P, init=0, sample_size=1000):
     sample_size : scalar(int), optional(default=1000)
         The length of the sample path.
 
+    random_state : scalar(int) or np.random.RandomState,
+                   optional(default=None)
+        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, int):
         X_0 = init
     else:
         cdf0 = np.cumsum(init)
-        u_0 = np.random.random(size=1)
+        u_0 = random_state.random_sample()
         X_0 = searchsorted(cdf0, u_0)
 
-    return MarkovChain(P).simulate(ts_length=sample_size, init=X_0)
+    mc = MarkovChain(P, random_state=random_state)
+    return mc.simulate(ts_length=sample_size, init=X_0)
diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index f8688c7fa..4755d2194 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -256,6 +256,27 @@ def test_simulate_for_matrices_with_C_F_orders():
     assert_array_equal(computed_F, sample_path)
 
 
+def test_mc_sample_path():
+    P = [[0.4, 0.6], [0.2, 0.8]]
+    init = 0
+    sample_size = 10
+
+    seed = 42
+
+    # init: integer
+    expected = [0, 0, 1, 1, 1, 0, 0, 0, 1, 1]
+    computed = mc_sample_path(P, init=init, sample_size=sample_size,
+                              random_state=seed)
+    assert_array_equal(computed, expected)
+
+    # init: distribution
+    distribution = (0.5, 0.5)
+    expected = [0, 1, 1, 1, 0, 0, 0, 1, 1, 1]
+    computed = mc_sample_path(P, init=distribution, sample_size=sample_size,
+                              random_state=seed)
+    assert_array_equal(computed, expected)
+
+
 @raises(ValueError)
 def test_raises_value_error_non_2dim():
     """Test with non 2dim input"""

From 0a752389f7154a311a17d0360e1668e12a208347 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 6 Jun 2015 20:46:44 +0900
Subject: [PATCH 05/12] MarkovChain.simulate refactored

---
 quantecon/markov/core.py | 58 ++++++++++++++++++----------------------
 1 file changed, 26 insertions(+), 32 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index a82528375..be4bd8139 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -313,12 +313,13 @@ def simulate(self, ts_length, init=None, num_reps=None):
         random_state = check_random_state(self.random_state)
 
         try:
-            num_reps = len(init)  # init is an array; make num_reps not None
+            k = len(init)  # init is an array
+            num_reps = k  # make num_reps not None
             init_states = np.asarray(init, dtype=int)
-        except:  # init is a scalar(int) or None
+        except TypeError:  # init is a scalar(int) or None
             k = 1 if num_reps is None else num_reps
             if init is None:
-                init_states = np.random.randint(self.n, size=k)
+                init_states = random_state.randint(self.n, size=k)
             elif isinstance(init, int):
                 init_states = np.ones(k, dtype=int) * init
             else:
@@ -326,8 +327,14 @@ def simulate(self, ts_length, init=None, num_reps=None):
                     'init must be int, array_like of ints, or None'
                 )
 
-        X = _simulate_markov_chain(self.cdfs, ts_length, init_states,
-                                   random_state=random_state)
+        # === set up array to store output === #
+        X = np.empty((k, ts_length), dtype=int)
+
+        # Random values, uniformly sampled from [0, 1)
+        random_values = random_state.random_sample(size=(k, ts_length-1))
+
+        # Generate sample paths and store in X
+        _generate_sample_paths(self.cdfs, init_states, random_values, out=X)
 
         if num_reps is None:
             return X[0]
@@ -335,54 +342,41 @@ def simulate(self, ts_length, init=None, num_reps=None):
             return X
 
 
-def _simulate_markov_chain(P_cdfs, ts_length, init_states, random_state):
+def _generate_sample_paths(P_cdfs, init_states, random_values, out):
     """
-    Main body of MarkovChain.simulate.
-
-    Generate len(init_states) sample paths of length ts_length.
+    Generate num_reps sample paths of length ts_length, where num_reps =
+    out.shape[0] and ts_length = out.shape[1].
 
     Parameters
     ----------
     P_cdfs : ndarray(float, ndim=2)
         Array containing as rows the CDFs of the state transition.
 
-    ts_length : scalar(int)
-        Length of each sample path.
-
     init_states : array_like(int, ndim=1)
-        Array containing the initial states.
+        Array containing the initial states. Its length must be equal to
+        num_reps.
 
-    random_state : np.random.RandomState
+    random_values : ndarray(float, ndim=2)
+        Array containing random values from [0, 1). Its shape must be
+        equal to (num_reps, ts_length-1)
 
-    Returns
-    -------
-    X : ndarray(int, ndim=2)
-        Array containing the sample paths, of shape (len(init_states),
-        ts_length)
+    out : ndarray(int, ndim=2)
+        Array to store the sample paths.
 
     Notes
     -----
     This routine is jit-complied if the module Numba is vailable.
 
     """
-    num_reps = len(init_states)
-
-    # === set up array to store output === #
-    X = np.empty((num_reps, ts_length), dtype=int)
-    X[:, 0] = init_states
+    num_reps, ts_length = out.shape
 
-    # Random values, uniformly sampled from [0, 1)
-    u = random_state.random_sample(size=(num_reps, ts_length-1))
-
-    # === generate the sample paths === #
     for i in range(num_reps):
+        out[i, 0] = init_states[i]
         for t in range(ts_length-1):
-            X[i, t+1] = searchsorted(P_cdfs[X[i, t]], u[i, t])
-
-    return X
+            out[i, t+1] = searchsorted(P_cdfs[out[i, t]], random_values[i, t])
 
 if numba_installed:
-    _simulate_markov_chain = jit(_simulate_markov_chain)
+    _generate_sample_paths = jit(nopython=True)(_generate_sample_paths)
 
 
 def mc_compute_stationary(P):

From 2f7b31cabe47d050d1689c22a3888cc00c66acf5 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 8 Aug 2015 13:45:51 +0900
Subject: [PATCH 06/12] Import statements corrected

---
 quantecon/markov/core.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index be4bd8139..21085aef1 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -2,7 +2,7 @@
 Authors: Chase Coleman, Spencer Lyon, Daisuke Oyama, Tom Sargent,
          John Stachurski
 
-Filename: mc_tools.py
+Filename: core.py
 
 This file contains some useful objects for handling a finite-state
 discrete-time Markov chain.
@@ -87,18 +87,23 @@
 from __future__ import division
 import numpy as np
 from fractions import gcd
-import sys
-from .graph_tools import DiGraph
 from .gth_solve import gth_solve
+from ..graph_tools import DiGraph
 
 # -Check if Numba is Available- #
 <<<<<<< HEAD
+<<<<<<< HEAD
 from ..util import searchsorted, numba_installed, jit
 =======
 from .external import numba_installed, jit
 
 from .util import searchsorted, check_random_state
 >>>>>>> `random_state` option, as well as tests, added
+=======
+from ..external import numba_installed, jit
+
+from ..util import searchsorted, check_random_state
+>>>>>>> Import statements corrected
 
 
 class MarkovChain(object):

From 83ebebfb686627b4d4869b82db71e10b69f7dbec Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 8 Aug 2015 14:36:27 +0900
Subject: [PATCH 07/12] `random_state` moved from `__init__` to `simulate`

based on discussion https://github.com/QuantEcon/QuantEcon.py/issues/146#issuecomment-116768263
---
 quantecon/markov/core.py | 28 +++++++++++++---------------
 1 file changed, 13 insertions(+), 15 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 21085aef1..8710ccfbc 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -118,13 +118,6 @@ class MarkovChain(object):
     P : array_like(float, ndim=2)
         The transition matrix.  Must be of shape n x n.
 
-    random_state : scalar(int) or np.random.RandomState,
-                   optional(default=None)
-        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.
-
     Attributes
     ----------
     P : see Parameters
@@ -160,7 +153,7 @@ class MarkovChain(object):
 
     """
 
-    def __init__(self, P, random_state=None):
+    def __init__(self, P):
         self.P = np.asarray(P)
 
         # Check Properties
@@ -179,9 +172,6 @@ def __init__(self, P, random_state=None):
         # The number of states
         self.n = self.P.shape[0]
 
-        # random state
-        self.random_state = random_state
-
         # To analyze the structure of P as a directed graph
         self._digraph = None
 
@@ -288,7 +278,7 @@ def cdfs(self):
             self._cdfs = cdfs
         return self._cdfs
 
-    def simulate(self, ts_length, init=None, num_reps=None):
+    def simulate(self, ts_length, init=None, num_reps=None, random_state=None):
         """
         Simulate time series of state transitions.
 
@@ -306,6 +296,13 @@ def simulate(self, ts_length, init=None, num_reps=None):
             Number of simulations. Relevant only when init is a scalar
             or None.
 
+        random_state : scalar(int) or np.random.RandomState,
+                       optional(default=None)
+            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 : ndarray(int, ndim=1 or 2)
@@ -315,7 +312,7 @@ def simulate(self, ts_length, init=None, num_reps=None):
             init is an array_like, otherwise k = num_reps.
 
         """
-        random_state = check_random_state(self.random_state)
+        random_state = check_random_state(random_state)
 
         try:
             k = len(init)  # init is an array
@@ -439,5 +436,6 @@ def mc_sample_path(P, init=0, sample_size=1000, random_state=None):
         u_0 = random_state.random_sample()
         X_0 = searchsorted(cdf0, u_0)
 
-    mc = MarkovChain(P, random_state=random_state)
-    return mc.simulate(ts_length=sample_size, init=X_0)
+    mc = MarkovChain(P)
+    return mc.simulate(ts_length=sample_size, init=X_0,
+                       random_state=random_state)

From d06a340726f4f0d2fb8e48fed3c6dcd92cbd8305 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 8 Aug 2015 16:15:10 +0900
Subject: [PATCH 08/12] test_simulate_lln added

---
 quantecon/markov/tests/test_core.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index 4755d2194..548ad29ef 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -11,7 +11,7 @@
 
 import numpy as np
 from numpy.testing import assert_allclose, assert_array_equal
-from nose.tools import eq_, raises
+from nose.tools import eq_, ok_, raises
 
 from quantecon.markov import (
     MarkovChain, mc_compute_stationary, mc_sample_path
@@ -233,6 +233,19 @@ def test_simulate_shape():
                            (num_reps, ts_length))
 
 
+def test_simulate_lln():
+    P = [[0.4, 0.6], [0.2, 0.8]]
+    stationary_dist = [0.25, 0.75]
+    mc = MarkovChain(P)
+
+    seed = 34
+    ts_length = 10**4
+    tol = 0.02
+
+    frequency_1 = mc.simulate(ts_length, random_state=seed).mean()
+    ok_(np.abs(frequency_1 - stationary_dist[1]) < tol)
+
+
 def test_simulate_for_matrices_with_C_F_orders():
     """
     Test MarkovChasin.simulate for matrices with C- and F-orders

From 80f2ffa61f33f6dadf8f86959b0ef1fad6d6144c Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 8 Aug 2015 17:04:39 +0900
Subject: [PATCH 09/12] test_simulate_ergodicity added

test_simulate_lln replaced by test_mc_sample_path_lln
---
 quantecon/markov/tests/test_core.py | 27 ++++++++++++++++++++++-----
 1 file changed, 22 insertions(+), 5 deletions(-)

diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index 548ad29ef..1347ba4a9 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -233,16 +233,19 @@ def test_simulate_shape():
                            (num_reps, ts_length))
 
 
-def test_simulate_lln():
+def test_simulate_ergodicity():
     P = [[0.4, 0.6], [0.2, 0.8]]
     stationary_dist = [0.25, 0.75]
+    init = 0
     mc = MarkovChain(P)
 
-    seed = 34
-    ts_length = 10**4
-    tol = 0.02
+    seed = 4433
+    ts_length = 100
+    num_reps = 300
+    tol = 0.1
 
-    frequency_1 = mc.simulate(ts_length, random_state=seed).mean()
+    x = mc.simulate(ts_length, init=init, num_reps=num_reps, random_state=seed)
+    frequency_1 = x[:, -1].mean()
     ok_(np.abs(frequency_1 - stationary_dist[1]) < tol)
 
 
@@ -290,6 +293,20 @@ def test_mc_sample_path():
     assert_array_equal(computed, expected)
 
 
+def test_mc_sample_path_lln():
+    P = [[0.4, 0.6], [0.2, 0.8]]
+    stationary_dist = [0.25, 0.75]
+    init = 0
+
+    seed = 4433
+    sample_size = 10**4
+    tol = 0.02
+
+    frequency_1 = mc_sample_path(P, init=init, sample_size=sample_size,
+                                 random_state=seed).mean()
+    ok_(np.abs(frequency_1 - stationary_dist[1]) < tol)
+
+
 @raises(ValueError)
 def test_raises_value_error_non_2dim():
     """Test with non 2dim input"""

From 7a6c227e31ec63279d81643419f608ecc9327e29 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Wed, 26 Aug 2015 21:45:59 +0900
Subject: [PATCH 10/12] Resolve conflicts

---
 quantecon/markov/core.py | 14 +-------------
 1 file changed, 1 insertion(+), 13 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 8710ccfbc..24813f363 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -91,19 +91,7 @@
 from ..graph_tools import DiGraph
 
 # -Check if Numba is Available- #
-<<<<<<< HEAD
-<<<<<<< HEAD
-from ..util import searchsorted, numba_installed, jit
-=======
-from .external import numba_installed, jit
-
-from .util import searchsorted, check_random_state
->>>>>>> `random_state` option, as well as tests, added
-=======
-from ..external import numba_installed, jit
-
-from ..util import searchsorted, check_random_state
->>>>>>> Import statements corrected
+from ..util import searchsorted, check_random_state, numba_installed, jit
 
 
 class MarkovChain(object):

From 45282cde8a687fba6987f624f4718b31eb1ae858 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Wed, 26 Aug 2015 23:03:11 +0900
Subject: [PATCH 11/12] simulate: Change the behavior for array init + num_reps

Now return len(init)*num_reps sample paths
---
 quantecon/markov/core.py            | 28 ++++++++++++++++++----------
 quantecon/markov/tests/test_core.py | 22 +++++++++++++++++++---
 2 files changed, 37 insertions(+), 13 deletions(-)

diff --git a/quantecon/markov/core.py b/quantecon/markov/core.py
index 24813f363..e3ee62d42 100644
--- a/quantecon/markov/core.py
+++ b/quantecon/markov/core.py
@@ -281,15 +281,14 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None):
             drawn.
 
         num_reps : scalar(int), optional(default=None)
-            Number of simulations. Relevant only when init is a scalar
-            or None.
+            Number of repetitions of simulation.
 
         random_state : scalar(int) or np.random.RandomState,
                        optional(default=None)
-            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.
+            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
         -------
@@ -297,17 +296,26 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None):
             Array containing the sample path(s), of shape (ts_length,)
             if init is a scalar (integer) or None and num_reps is None;
             of shape (k, ts_length) otherwise, where k = len(init) if
-            init is an array_like, otherwise k = num_reps.
+            (init, num_reps) = (array, None), k = num_reps if (init,
+            num_reps) = (int or None, int), and k = len(init)*num_reps
+            if (init, num_reps) = (array, int).
 
         """
         random_state = check_random_state(random_state)
+        dim = 1  # Dimension of the returned array: 1 or 2
 
         try:
             k = len(init)  # init is an array
-            num_reps = k  # make num_reps not None
+            dim = 2
             init_states = np.asarray(init, dtype=int)
+            if num_reps is not None:
+                k *= num_reps
+                init_states = np.tile(init_states, num_reps)
         except TypeError:  # init is a scalar(int) or None
-            k = 1 if num_reps is None else num_reps
+            k = 1
+            if num_reps is not None:
+                dim = 2
+                k = num_reps
             if init is None:
                 init_states = random_state.randint(self.n, size=k)
             elif isinstance(init, int):
@@ -326,7 +334,7 @@ def simulate(self, ts_length, init=None, num_reps=None, random_state=None):
         # Generate sample paths and store in X
         _generate_sample_paths(self.cdfs, init_states, random_values, out=X)
 
-        if num_reps is None:
+        if dim == 1:
             return X[0]
         else:
             return X
diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index 1347ba4a9..1d482cb44 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -224,15 +224,31 @@ def test_simulate_shape():
     assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
                        (ts_length,))
 
-    for (ts_length, init, num_reps) in [(10, [0, 1], None), (10, [0, 1], 3)]:
-        assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
-                           (len(init), ts_length))
+    (ts_length, init, num_reps) = (10, [0, 1], None)
+    assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
+                       (len(init), ts_length))
+
+    (ts_length, init, num_reps) = (10, [0, 1], 3)
+    assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
+                       (len(init)*num_reps, ts_length))
 
     for (ts_length, init, num_reps) in [(10, None, 3), (10, None, 1)]:
         assert_array_equal(mc.simulate(ts_length, init, num_reps).shape,
                            (num_reps, ts_length))
 
 
+def test_simulate_init_array_num_reps():
+    P = [[0.4, 0.6], [0.2, 0.8]]
+    mc = MarkovChain(P)
+
+    ts_length = 10
+    init=[0, 1]
+    num_reps=3
+
+    X = mc.simulate(ts_length, init, num_reps)
+    assert_array_equal(X[:, 0], init*num_reps)
+
+
 def test_simulate_ergodicity():
     P = [[0.4, 0.6], [0.2, 0.8]]
     stationary_dist = [0.25, 0.75]

From c535ab7ca089700ce525ef0f4e0146dcdd5159d6 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Wed, 26 Aug 2015 23:03:51 +0900
Subject: [PATCH 12/12] Remove unused import

---
 quantecon/markov/tests/test_core.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/quantecon/markov/tests/test_core.py b/quantecon/markov/tests/test_core.py
index 1d482cb44..736e42f1a 100644
--- a/quantecon/markov/tests/test_core.py
+++ b/quantecon/markov/tests/test_core.py
@@ -17,8 +17,6 @@
     MarkovChain, mc_compute_stationary, mc_sample_path
 )
 
-from quantecon.util import numba_installed, jit
-
 
 def list_of_array_equal(s, t):
     """