Repository: econ-ark/HARK

Issue:
discrete distributions should get a random seed on initialization
`DiscreteDistribution.drawDiscrete` takes a `seed` argument.

Instead, this seed can be an argument to the DiscreteDistribution upon initialization.

Discussion:

Patch:
diff --git a/HARK/ConsumptionSaving/ConsAggShockModel.py b/HARK/ConsumptionSaving/ConsAggShockModel.py
index 1fb54229..7b04d578 100644
--- a/HARK/ConsumptionSaving/ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/ConsAggShockModel.py
@@ -480,8 +480,7 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
                 
                 # Get random draws of income shocks from the discrete distribution
                 ShockDraws = IncomeDstnNow.drawDiscrete(N,
-                                                        exact_match=True,
-                                                        seed=self.RNG.randint(0, 2**31-1))
+                                                        exact_match=True)
                 # Permanent "shock" includes expected growth
                 PermShkNow[these] = ShockDraws[0]*PermGroFacNow
                 TranShkNow[these] = ShockDraws[1]
@@ -496,8 +495,7 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
             
             # Get random draws of income shocks from the discrete distribution
             ShockDraws = IncomeDstnNow.drawDiscrete(N,
-                                                    exact_match=True,
-                                                    seed=self.RNG.randint(0, 2**31-1))
+                                                    exact_match=True)
             # Permanent "shock" includes expected growth
             PermShkNow[these] = ShockDraws[0]*PermGroFacNow
             TranShkNow[these] = ShockDraws[1]
@@ -1133,8 +1131,7 @@ class CobbDouglasEconomy(Market):
         Events = np.arange(self.AggShkDstn.pmf.size)  # just a list of integers
         EventDraws = self.AggShkDstn.drawDiscrete(
             N=sim_periods,
-            X=Events,
-            seed=0)
+            X=Events)
         PermShkAggHist = self.AggShkDstn.X[0][EventDraws]
         TranShkAggHist = self.AggShkDstn.X[1][EventDraws]
 
@@ -1369,8 +1366,7 @@ class SmallOpenEconomy(Market):
         Events = np.arange(self.AggShkDstn.pmf.size)  # just a list of integers
         EventDraws = self.AggShkDstn.drawDiscrete(
             N=sim_periods,
-            X=Events,
-            seed=0)
+            X=Events)
         PermShkAggHist = self.AggShkDstn.X[0][EventDraws]
         TranShkAggHist = self.AggShkDstn.X[1][EventDraws]
 
@@ -1561,7 +1557,7 @@ class CobbDouglasMarkovEconomy(CobbDouglasEconomy):
         PermShkAggHistAll = np.zeros((StateCount, sim_periods))
         TranShkAggHistAll = np.zeros((StateCount, sim_periods))
         for i in range(StateCount):
-            AggShockDraws = self.AggShkDstn[i].drawDiscrete(N=sim_periods, seed=0)
+            AggShockDraws = self.AggShkDstn[i].drawDiscrete(N=sim_periods)
             PermShkAggHistAll[i, :] = AggShockDraws[0,:]
             TranShkAggHistAll[i, :] = AggShockDraws[1,:]
 
@@ -1628,7 +1624,7 @@ class CobbDouglasMarkovEconomy(CobbDouglasEconomy):
 
         # Add histories until each state has been visited at least state_T_min times
         while go:
-            draws = Uniform().draw(N=self.act_T_orig, seed=loops)
+            draws = Uniform(seed=loops).draw(N=self.act_T_orig)
             for s in range(draws.size):  # Add act_T_orig more periods
                 MrkvNow_hist[t] = MrkvNow
                 MrkvNow = np.searchsorted(cutoffs[MrkvNow, :], draws[s])
diff --git a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
index ed4bf448..f40c394f 100644
--- a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
@@ -1054,14 +1054,15 @@ class GenIncProcessConsumerType(IndShockConsumerType):
         # Simulate the distribution of persistent income levels by t_cycle in a lifecycle model
         if self.cycles == 1:
             pLvlNow = Lognormal(self.pLvlInitMean,
-                                sigma=self.pLvlInitStd,).draw(self.AgentCount, seed=31382)
+                                sigma=self.pLvlInitStd,
+                                seed=31382
+            ).draw(self.AgentCount)
             pLvlGrid = []  # empty list of time-varying persistent income grids
             # Calculate distribution of persistent income in each period of lifecycle
             for t in range(len(self.PermShkStd)):
                 if t > 0:
                     PermShkNow = self.PermShkDstn[t-1].drawDiscrete(
-                        N=self.AgentCount,
-                        seed=t)
+                        N=self.AgentCount)
                     pLvlNow = self.pLvlNextFunc[t-1](pLvlNow)*PermShkNow
                 pLvlGrid.append(getPercentiles(pLvlNow, percentiles=self.pLvlPctiles))
 
@@ -1069,21 +1070,22 @@ class GenIncProcessConsumerType(IndShockConsumerType):
         elif self.cycles == 0:
             T_long = 1000  # Number of periods to simulate to get to "stationary" distribution
             pLvlNow = Lognormal(mu=self.pLvlInitMean,
-                                sigma=self.pLvlInitStd).draw(self.AgentCount, seed=31382)
+                                sigma=self.pLvlInitStd,
+                                seed=31382).draw(self.AgentCount)
             t_cycle = np.zeros(self.AgentCount, dtype=int)
             for t in range(T_long):
                 LivPrb = LivPrbAll[t_cycle]  # Determine who dies and replace them with newborns
-                draws = Uniform().draw(self.AgentCount, seed=t)
+                draws = Uniform(seed=t).draw(self.AgentCount)
                 who_dies = draws > LivPrb
                 pLvlNow[who_dies] = Lognormal(self.pLvlInitMean,
-                                              self.pLvlInitStd).draw(np.sum(who_dies),  seed=t+92615)
+                                              self.pLvlInitStd,
+                                              seed=t+92615).draw(np.sum(who_dies))
                 t_cycle[who_dies] = 0
                 
                 for j in range(self.T_cycle):  # Update persistent income
                     these = t_cycle == j
                     PermShkTemp = self.PermShkDstn[j].drawDiscrete(
-                        N=np.sum(these),
-                        seed=t+13*j)
+                        N=np.sum(these))
                     pLvlNow[these] = self.pLvlNextFunc[j](pLvlNow[these])*PermShkTemp
                 t_cycle = t_cycle + 1
                 t_cycle[t_cycle == self.T_cycle] = 0
@@ -1120,12 +1122,12 @@ class GenIncProcessConsumerType(IndShockConsumerType):
         # Get and store states for newly born agents
         N = np.sum(which_agents)  # Number of new consumers to make
         aNrmNow_new = Lognormal(self.aNrmInitMean,
-                                self.aNrmInitStd).draw(
-                                    N,
-                                    seed=self.RNG.randint(0, 2**31-1))
+                                self.aNrmInitStd,
+                                seed=self.RNG.randint(0, 2**31-1)).draw(
+                                    N)
         self.pLvlNow[which_agents] = Lognormal(self.pLvlInitMean,
-                                               self.pLvlInitStd).draw(N, 
-                                                                       seed=self.RNG.randint(0, 2**31-1))
+                                               self.pLvlInitStd, 
+                                               seed=self.RNG.randint(0, 2**31-1)).draw(N)
         self.aLvlNow[which_agents] = aNrmNow_new*self.pLvlNow[which_agents]
         self.t_age[which_agents] = 0  # How many periods since each agent was born
         self.t_cycle[which_agents] = 0  # Which period of the cycle each agent is currently in
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index 2bda5cf3..0488a846 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -1574,11 +1574,15 @@ class PerfForesightConsumerType(AgentType):
         '''
         # Get and store states for newly born agents
         N = np.sum(which_agents) # Number of new consumers to make
-        self.aNrmNow[which_agents] = Lognormal(mu=self.aNrmInitMean,
-                                               sigma=self.aNrmInitStd).draw(N,seed=self.RNG.randint(0,2**31-1))
+        self.aNrmNow[which_agents] = Lognormal(
+            mu=self.aNrmInitMean,
+            sigma=self.aNrmInitStd,
+            seed=self.RNG.randint(0,2**31-1)).draw(N)
         pLvlInitMeanNow = self.pLvlInitMean + np.log(self.PlvlAggNow) # Account for newer cohorts having higher permanent income
-        self.pLvlNow[which_agents] = Lognormal(pLvlInitMeanNow,
-                                               self.pLvlInitStd,).draw(N,seed=self.RNG.randint(0,2**31-1))
+        self.pLvlNow[which_agents] = Lognormal(
+            pLvlInitMeanNow,
+            self.pLvlInitStd,
+            seed=self.RNG.randint(0,2**31-1)).draw(N)
         self.t_age[which_agents]   = 0 # How many periods since each agent was born
         self.t_cycle[which_agents] = 0 # Which period of the cycle each agent is currently in
         return None
@@ -1600,7 +1604,8 @@ class PerfForesightConsumerType(AgentType):
         # Determine who dies
         DiePrb_by_t_cycle = 1.0 - np.asarray(self.LivPrb)
         DiePrb = DiePrb_by_t_cycle[self.t_cycle-1] # Time has already advanced, so look back one
-        DeathShks = Uniform().draw(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        DeathShks = Uniform(
+            seed=self.RNG.randint(0,2**31-1)).draw(N=self.AgentCount)
         which_agents = DeathShks < DiePrb
         if self.T_age is not None: # Kill agents that have lived for too many periods
             too_old = self.t_age >= self.T_age
@@ -2019,8 +2024,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
                 IncomeDstnNow    = self.IncomeDstn[t-1] # set current income distribution
                 PermGroFacNow    = self.PermGroFac[t-1] # and permanent growth factor
                 # Get random draws of income shocks from the discrete distribution
-                EventDraws       = IncomeDstnNow.draw_events(N,
-                                                             seed=self.RNG.randint(0,2**31-1))
+                EventDraws       = IncomeDstnNow.draw_events(N)
 
                 PermShkNow[these] = IncomeDstnNow.X[0][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
                 TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
@@ -2034,8 +2038,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
             PermGroFacNow    = self.PermGroFac[0] # and permanent growth factor
            
             # Get random draws of income shocks from the discrete distribution
-            EventDraws       = IncomeDstnNow.draw_events(N,
-                                                         seed=self.RNG.randint(0,2**31-1))
+            EventDraws       = IncomeDstnNow.draw_events(N)
             PermShkNow[these] = IncomeDstnNow.X[0][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
             TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
 #        PermShkNow[newborn] = 1.0
diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py b/HARK/ConsumptionSaving/ConsMarkovModel.py
index d7386f69..55469b09 100644
--- a/HARK/ConsumptionSaving/ConsMarkovModel.py
+++ b/HARK/ConsumptionSaving/ConsMarkovModel.py
@@ -778,7 +778,7 @@ class MarkovConsumerType(IndShockConsumerType):
     def initializeSim(self):
         IndShockConsumerType.initializeSim(self)
         if self.global_markov:  #Need to initialize markov state to be the same for all agents
-            base_draw = Uniform().draw(1,seed=self.RNG.randint(0,2**31-1))
+            base_draw = Uniform(seed=self.RNG.randint(0,2**31-1)).draw(1)
             Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit))
             self.MrkvNow = np.ones(self.AgentCount)*np.searchsorted(Cutoffs,base_draw).astype(int)
         self.MrkvNow = self.MrkvNow.astype(int)
@@ -800,7 +800,7 @@ class MarkovConsumerType(IndShockConsumerType):
         # Determine who dies
         LivPrb = np.array(self.LivPrb)[self.t_cycle-1,self.MrkvNow] # Time has already advanced, so look back one
         DiePrb = 1.0 - LivPrb
-        DeathShks = Uniform().draw(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        DeathShks = Uniform(seed=self.RNG.randint(0,2**31-1)).draw(N=self.AgentCount)
         which_agents = DeathShks < DiePrb
         if self.T_age is not None: # Kill agents that have lived for too many periods
             too_old = self.t_age >= self.T_age
@@ -824,7 +824,7 @@ class MarkovConsumerType(IndShockConsumerType):
         IndShockConsumerType.simBirth(self,which_agents) # Get initial assets and permanent income
         if not self.global_markov:  #Markov state is not changed if it is set at the global level
             N = np.sum(which_agents)
-            base_draws = Uniform().draw(N,seed=self.RNG.randint(0,2**31-1))
+            base_draws = Uniform(seed=self.RNG.randint(0,2**31-1)).draw(N)
             Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit))
             self.MrkvNow[which_agents] = np.searchsorted(Cutoffs,base_draws).astype(int)
             
@@ -844,9 +844,11 @@ class MarkovConsumerType(IndShockConsumerType):
         '''
         # Draw random numbers that will be used to determine the next Markov state
         if self.global_markov:
-            base_draws = np.ones(self.AgentCount)*Uniform().draw(1,seed=self.RNG.randint(0,2**31-1))
+            base_draws = np.ones(self.AgentCount)*Uniform(
+                seed=self.RNG.randint(0,2**31-1)).draw(1)
         else:
-            base_draws = Uniform().draw(self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+            base_draws = Uniform(
+                seed=self.RNG.randint(0,2**31-1)).draw(self.AgentCount)
         dont_change = self.t_age == 0 # Don't change Markov state for those who were just born (unless global_markov)
         if self.t_sim == 0: # Respect initial distribution of Markov states
             dont_change[:] = True
@@ -900,8 +902,7 @@ class MarkovConsumerType(IndShockConsumerType):
 
                     # Get random draws of income shocks from the discrete distribution
                     EventDraws       = IncomeDstnNow.draw_events(
-                        N,
-                        seed=self.RNG.randint(0,2**31-1))
+                        N)
                     PermShkNow[these] = IncomeDstnNow.X[0][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
                     TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
         newborn = self.t_age == 0
diff --git a/HARK/ConsumptionSaving/ConsRepAgentModel.py b/HARK/ConsumptionSaving/ConsRepAgentModel.py
index 577a7c84..754bdc2a 100644
--- a/HARK/ConsumptionSaving/ConsRepAgentModel.py
+++ b/HARK/ConsumptionSaving/ConsRepAgentModel.py
@@ -304,7 +304,7 @@ class RepAgentMarkovConsumerType(RepAgentConsumerType):
         None
         '''
         cutoffs = np.cumsum(self.MrkvArray[self.MrkvNow,:])
-        MrkvDraw = Uniform().draw(N=1,seed=self.RNG.randint(0,2**31-1))
+        MrkvDraw = Uniform(seed=self.RNG.randint(0,2**31-1)).draw(N=1)
         self.MrkvNow = np.searchsorted(cutoffs,MrkvDraw)
 
         t = self.t_cycle[0]
@@ -312,8 +312,7 @@ class RepAgentMarkovConsumerType(RepAgentConsumerType):
         IncomeDstnNow    = self.IncomeDstn[t-1][i] # set current income distribution
         PermGroFacNow    = self.PermGroFac[t-1][i] # and permanent growth factor
         # Get random draws of income shocks from the discrete distribution
-        EventDraw        =         IncomeDstnNow.draw_events(1,
-                                                             seed=self.RNG.randint(0,2**31-1))
+        EventDraw        =         IncomeDstnNow.draw_events(1)
         PermShkNow = IncomeDstnNow.X[0][EventDraw]*PermGroFacNow # permanent "shock" includes expected growth
         TranShkNow = IncomeDstnNow.X[1][EventDraw]
         self.PermShkNow = np.array(PermShkNow)
diff --git a/HARK/distribution.py b/HARK/distribution.py
index 4e862785..542ac559 100644
--- a/HARK/distribution.py
+++ b/HARK/distribution.py
@@ -5,16 +5,40 @@ from scipy.special import erf, erfc
 import scipy.stats as stats
 import warnings
 
+
+class Distribution():
+
+    def __init__(self, seed = 0):
+        '''
+        Initialize the distribution.
+
+        Parameters
+        ----------
+        seed : int
+            Seed for random number generator.
+        '''
+        self.RNG = np.random.RandomState(seed)
+        self.seed = seed
+
+    def reset(self):
+        '''
+        Reset the random number generator of this distribution.
+
+        Parameters
+        ----------
+        '''
+        self.RNG = np.random.RandomState(self.seed)
+    
 ### CONTINUOUS DISTRIBUTIONS
 
-class Lognormal():
+class Lognormal(Distribution):
     """
     A Lognormal distribution
     """
     mu = None
     sigma = None
 
-    def __init__(self, mu = 0.0, sigma = 1.0):
+    def __init__(self, mu = 0.0, sigma = 1.0, seed = 0):
         '''
         Initialize the distribution.
 
@@ -26,11 +50,15 @@ class Lognormal():
         sigma : float or [float]
             One or more standard deviations. Number of elements T in sigma
             determines number of rows of output.
+        seed : int
+            Seed for random number generator.
         '''
         self.mu = mu
         self.sigma = sigma
+        # Set up the RNG
+        super().__init__(seed)
 
-    def draw(self, N, seed=0):
+    def draw(self, N):
         '''
         Generate arrays of lognormal draws. The sigma input can be a number
         or list-like.  If a number, output is a length N array of draws from the
@@ -42,8 +70,6 @@ class Lognormal():
         ----------
         N : int
             Number of draws in each row.
-        seed : int
-            Seed for random number generator.
 
         Returns:
         ------------
@@ -51,25 +77,23 @@ class Lognormal():
             T-length list of arrays of mean one lognormal draws each of size N, or
             a single array of size N (if sigma is a scalar).
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
 
         if isinstance(self.sigma,float): # Return a single array of length N
             if self.sigma == 0:
                 draws = np.exp(self.mu)*np.ones(N)
             else:
-                draws = RNG.lognormal(mean=self.mu,
-                                      sigma=self.sigma,
-                                      size=N)
+                draws = self.RNG.lognormal(mean=self.mu,
+                                           sigma=self.sigma,
+                                           size=N)
         else: # Set up empty list to populate, then loop and populate list with draws
             draws=[]
             for j in range(len(self.sigma)):
                 if self.sigma[j] == 0:
                     draws.append(np.exp(self.mu[j])*np.ones(N))
                 else:
-                    draws.append(RNG.lognormal(mean=self.mu[j],
-                                               sigma=self.sigma[j],
-                                               size=N))
+                    draws.append(self.RNG.lognormal(mean=self.mu[j],
+                                                    sigma=self.sigma[j],
+                                                    size=N))
         return draws
 
     def approx(self, N, tail_N=0, tail_bound=None, tail_order=np.e):
@@ -148,23 +172,27 @@ class Lognormal():
         else:
             pmf = np.ones(N)/N
             X   = np.exp(self.mu)*np.ones(N)
-        return DiscreteDistribution(pmf, X)
+        return DiscreteDistribution(pmf,
+                                    X,
+                                    seed=self.RNG.randint(0,
+                                                          2**31 - 1,
+                                                          dtype='int32'))
 
 class MeanOneLogNormal(Lognormal):
 
-    def __init__(self, sigma=1.0):
+    def __init__(self, sigma=1.0, seed=0):
         mu = -0.5*sigma**2
-        super().__init__(mu=mu, sigma=sigma)
+        super().__init__(mu=mu, sigma=sigma, seed=seed)
 
 
-class Normal():
+class Normal(Distribution):
     """
     A Normal distribution.
     """
     mu = None
     sigma = None
 
-    def __init__(self, mu = 0.0, sigma = 1.0):
+    def __init__(self, mu = 0.0, sigma = 1.0, seed=0):
         '''
         Initialize the distribution.
 
@@ -176,9 +204,12 @@ class Normal():
         sigma : float or [float]
             One or more standard deviations. Number of elements T in sigma
             determines number of rows of output.
+        seed : int
+            Seed for random number generator.
         '''
         self.mu = mu
         self.sigma = sigma
+        super().__init__(seed)
 
     def draw(self, N, seed=0):
         '''
@@ -192,8 +223,6 @@ class Normal():
         ----------
         N : int
             Number of draws in each row.
-        seed : int
-            Seed for random number generator.
 
         Returns
         -------
@@ -201,15 +230,12 @@ class Normal():
             T-length list of arrays of normal draws each of size N, or a single array
             of size N (if sigma is a scalar).
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
-
         if isinstance(self.sigma,float): # Return a single array of length N
-            draws = self.sigma*RNG.randn(N) + self.mu
+            draws = self.sigma*self.RNG.randn(N) + self.mu
         else: # Set up empty list to populate, then loop and populate list with draws
             draws=[]
             for t in range(len(sigma)):
-                draws.append(self.sigma[t]*RNG.randn(N) + self.mu[t])
+                draws.append(self.sigma[t]*self.RNG.randn(N) + self.mu[t])
         return draws
 
     def approx(self, N):
@@ -221,9 +247,13 @@ class Normal():
         pmf = w*np.pi**-0.5
         # correct x
         X = math.sqrt(2.0)*self.sigma*x + self.mu
-        return DiscreteDistribution(pmf, X)
+        return DiscreteDistribution(pmf,
+                                    X,
+                                    seed=self.RNG.randint(0,
+                                                          2**31 - 1,
+                                                          dtype='int32'))
 
-class Weibull():
+class Weibull(Distribution):
     '''
     A Weibull distribution
     '''
@@ -231,23 +261,28 @@ class Weibull():
     scale = None
     shape = None
 
-    def __init__(self, scale=1.0, shape=1.0):
+    def __init__(self, scale=1.0, shape=1.0, seed=0):
         '''
         Initialize the distribution.
 
         Parameters
         ----------
         scale : float or [float]
-            One or more scales.  Number of elements T in scale determines number of
-        rows of output.
+            One or more scales.  Number of elements T in scale
+            determines number of
+            rows of output.
         shape : float or [float]
             One or more shape parameters. Number of elements T in scale
             determines number of rows of output.
+        seed : int
+            Seed for random number generator.
         '''
         self.scale = scale
         self.shape = shape
+        # Set up the RNG
+        super().__init__(seed)
 
-    def draw(self, N, seed=0):
+    def draw(self, N):
         '''
         Generate arrays of Weibull draws.  The scale and shape inputs can be
         numbers or list-likes.  If a number, output is a length N array of draws from
@@ -263,8 +298,6 @@ class Weibull():
         ----------
         N : int
             Number of draws in each row.
-        seed : int
-            Seed for random number generator.
 
         Returns:
         ------------
@@ -272,20 +305,17 @@ class Weibull():
             T-length list of arrays of Weibull draws each of size N, or a single
             array of size N (if sigma is a scalar).
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
-
         if self.scale == 1:
             scale = float(self.scale)
         if isinstance(self.scale,float): # Return a single array of length N
-            draws = self.scale*(-np.log(1.0-RNG.rand(N)))**(1.0/self.shape)
+            draws = self.scale*(-np.log(1.0-self.RNG.rand(N)))**(1.0/self.shape)
         else: # Set up empty list to populate, then loop and populate list with draws
             draws=[]
             for t in range(len(self.scale)):
-                draws.append(self.scale[t]*(-np.log(1.0-RNG.rand(N)))**(1.0/self.shape[t]))
+                draws.append(self.scale[t]*(-np.log(1.0-self.RNG.rand(N)))**(1.0/self.shape[t]))
         return draws
 
-class Uniform():
+class Uniform(Distribution):
     """
     A Uniform distribution.
     """
@@ -293,23 +323,29 @@ class Uniform():
     bot = None
     top = None
 
-    def __init__(self, bot = 0.0, top = 1.0):
+    def __init__(self, bot = 0.0, top = 1.0, seed=0):
         '''
         Initialize the distribution.
 
         Parameters
         ----------
         bot : float or [float]
-            One or more bottom values.  Number of elements T in mu determines number
+            One or more bottom values.
+            Number of elements T in mu determines number
             of rows of output.
         top : float or [float]
-            One or more top values. Number of elements T in top determines number of
+            One or more top values.
+            Number of elements T in top determines number of
             rows of output.
+        seed : int
+            Seed for random number generator.
         '''
         self.bot = bot
         self.top = top
+        # Set up the RNG
+        self.RNG = np.random.RandomState(seed)
 
-    def draw(self, N, seed=0):
+    def draw(self, N):
         '''
         Generate arrays of uniform draws.  The bot and top inputs can be numbers or
         list-likes.  If a number, output is a length N array of draws from the
@@ -321,8 +357,6 @@ class Uniform():
         ----------
         N : int
             Number of draws in each row.
-        seed : int
-            Seed for random number generator.
 
         Returns
         -------
@@ -330,15 +364,12 @@ class Uniform():
             T-length list of arrays of uniform draws each of size N, or a single
             array of size N (if sigma is a scalar).
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
-
         if isinstance(self.bot,float) or isinstance(self.bot,int): # Return a single array of size N
-            draws = self.bot + (self.top - self.bot)*RNG.rand(N)
+            draws = self.bot + (self.top - self.bot)*self.RNG.rand(N)
         else: # Set up empty list to populate, then loop and populate list with draws
             draws=[]
             for t in range(len(bot)):
-                draws.append(self.bot[t] + (self.top[t] - self.bot[t])*RNG.rand(N))
+                draws.append(self.bot[t] + (self.top[t] - self.bot[t])*self.RNG.rand(N))
         return draws
 
     def approx(self, N):
@@ -360,19 +391,24 @@ class Uniform():
         center = (self.top+self.bot)/2.0
         width = (self.top-self.bot)/2.0
         X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0)
-        return DiscreteDistribution(pmf,X)
+        return DiscreteDistribution(pmf,
+                                    X,
+                                    seed=self.RNG.randint(0,
+                                                          2**31 - 1,
+                                                          dtype='int32')
+                                    )
 
 
 ### DISCRETE DISTRIBUTIONS
 
-class Bernoulli():
+class Bernoulli(Distribution):
     """
     A Bernoulli distribution.
     """
 
     p = None
 
-    def __init__(self, p = 0.5):
+    def __init__(self, p = 0.5, seed = 0):
         '''
         Initialize the distribution.
 
@@ -380,10 +416,16 @@ class Bernoulli():
         ----------
         p : float or [float]
             Probability or probabilities of the event occurring (True).
+
+        seed : int
+            Seed for random number generator.
         '''
         self.p = p
+        # Set up the RNG
+        super().__init__(seed)
+
 
-    def draw(self, N, seed = 0):
+    def draw(self, N):
         '''
         Generates arrays of booleans drawn from a simple Bernoulli distribution.
         The input p can be a float or a list-like of floats; its length T determines
@@ -394,8 +436,6 @@ class Bernoulli():
         ---------
         N : int
             Number of draws in each row.
-        seed : int
-            Seed for random number generator.
 
         Returns
         -------
@@ -403,20 +443,15 @@ class Bernoulli():
             T-length list of arrays of Bernoulli draws each of size N, or a single
         array of size N (if sigma is a scalar).
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
-
         if isinstance(self.p,float):# Return a single array of size N
-            draws = RNG.uniform(size=N) < self.p
+            draws = self.RNG.uniform(size=N) < self.p
         else: # Set up empty list to populate, then loop and populate list with draws:
             draws=[]
             for t in range(len(self.p)):
-                draws.append(RNG.uniform(size=N) < self.p[t])
+                draws.append(self.RNG.uniform(size=N) < self.p[t])
         return draws
 
-
-
-class DiscreteDistribution():
+class DiscreteDistribution(Distribution):
     """
     A representation of a discrete probability distribution.
 
@@ -424,7 +459,7 @@ class DiscreteDistribution():
     pmf = None
     X = None
 
-    def __init__(self, pmf, X):
+    def __init__(self, pmf, X, seed=0):
         '''
         Initialize a discrete distribution.
 
@@ -435,10 +470,14 @@ class DiscreteDistribution():
         X : np.array or [np.array]
             Discrete point values for each probability mass.
             May be multivariate (list of arrays).
-        
+        seed : int
+            Seed for random number generator.
         '''
         self.pmf = pmf
         self.X = X
+        # Set up the RNG
+        super().__init__(seed)
+        
 
         # Very quick and incomplete parameter check:
         # TODO: Check that pmf and X arrays have same length.
@@ -449,15 +488,13 @@ class DiscreteDistribution():
         else:
             return 1
 
-    def draw_events(self, n, seed=0):
+    def draw_events(self, n):
         '''
         Draws N 'events' from the distribution PMF.
         These events are indices into X.
         '''
-        # Set up the RNG
-        RNG = np.random.RandomState(seed)
         # Generate a cumulative distribution
-        base_draws = RNG.uniform(size=n)
+        base_draws = self.RNG.uniform(size=n)
         cum_dist = np.cumsum(self.pmf)
 
         # Convert the basic uniform draws into discrete draws
@@ -465,7 +502,7 @@ class DiscreteDistribution():
 
         return indices
 
-    def drawDiscrete(self, N,X=None,exact_match=False,seed=0):
+    def drawDiscrete(self, N,X=None,exact_match=False):
         '''
         Simulates N draws from a discrete distribution with probabilities P and outcomes X.
 
@@ -483,8 +520,6 @@ class DiscreteDistribution():
             a random permutation of the N-length list that best fits the discrete
             distribution.  When False (default), each draw is independent from the
             others and the result could deviate from the input.
-        seed : int
-            Seed for random number generator.
 
         Returns
         -------
@@ -501,15 +536,10 @@ class DiscreteDistribution():
             X = X
             J = 1
 
-        # Draw indices whose empirical distribution closely matches discrete pmf
         if exact_match:
-            # Set up the RNG
-            RNG = np.random.RandomState(seed)
-            
             events = np.arange(self.pmf.size) # just a list of integers
             cutoffs = np.round(np.cumsum(self.pmf)*N).astype(int) # cutoff points between discrete outcomes
             top = 0
-            
             # Make a list of event indices that closely matches the discrete distribution
             event_list        = []
             for j in range(events.size):
@@ -518,11 +548,11 @@ class DiscreteDistribution():
                 event_list += (top-bot)*[events[j]]
 
             # Randomly permute the event indices
-            indices = RNG.permutation(event_list)
+            indices = self.RNG.permutation(event_list)
         
         # Draw event indices randomly from the discrete distribution
         else:
-            indices = self.draw_events(N, seed=seed)
+            indices = self.draw_events(N)
             
         # Create and fill in the output array of draws based on the output of event indices
         if J > 1:
@@ -534,9 +564,9 @@ class DiscreteDistribution():
 
         return draws
 
-def approxLognormalGaussHermite(N, mu=0.0, sigma=1.0):
+def approxLognormalGaussHermite(N, mu=0.0, sigma=1.0, seed=0):
     d = Normal(mu, sigma).approx(N)
-    return DiscreteDistribution(d.pmf, np.exp(d.X))
+    return DiscreteDistribution(d.pmf, np.exp(d.X), seed=seed)
 
 
 def calcNormalStyleParsFromLognormalPars(avgLognormal, stdLognormal):


Test patch:
diff --git a/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py b/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
index 43919f31..7a17a3c4 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
@@ -45,7 +45,7 @@ class testAggShockConsumerType(unittest.TestCase):
         
         self.economy.AFunc = self.economy.dynamics.AFunc
         self.assertAlmostEqual(self.economy.AFunc.slope,
-                         1.124330884813638)
+                               1.1374781287418916)
 
 
 class testAggShockMarkovConsumerType(unittest.TestCase):
@@ -79,5 +79,5 @@ class testAggShockMarkovConsumerType(unittest.TestCase):
         
         self.economy.AFunc = self.economy.dynamics.AFunc
         self.assertAlmostEqual(self.economy.AFunc[0].slope,
-                               1.0801777346256896)
+                               1.0897651487240472)
                          
diff --git a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
index 5a98c3ef..374803b1 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
@@ -93,4 +93,4 @@ class testPersistentShockConsumerType(unittest.TestCase):
 
         self.assertAlmostEqual(
             np.mean(self.agent.history['mLvlNow']),
-            1.2043701902476343)
+            1.2063544810887799)
diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
index 59f1bad2..9353cba4 100644
--- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
+++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
@@ -23,11 +23,11 @@ class testIndShockConsumerType(unittest.TestCase):
         self.agent.getShocks()
 
         self.assertEqual(self.agent.PermShkNow[0],
-                         1.0050166461586711)
+                         0.9686755529883603)
         self.assertEqual(self.agent.PermShkNow[1],
                          1.0050166461586711)
         self.assertEqual(self.agent.TranShkNow[0],
-                         1.1176912196531754)
+                         1.2093790234554653)
 
     def test_ConsIndShockSolverBasic(self):
         LifecycleExample = IndShockConsumerType(
@@ -105,10 +105,10 @@ class testIndShockConsumerType(unittest.TestCase):
         print(self.agent.aLvlNow)
 
         self.assertAlmostEqual(self.agent.MPCnow[1],
-                               0.5535801655448935)
+                               0.5711503906043797)
 
         self.assertAlmostEqual(self.agent.aLvlNow[1],
-                               0.18832361)
+                               0.21251164208206255)
 
 
 class testBufferStock(unittest.TestCase):
diff --git a/HARK/tests/test_distribution.py b/HARK/tests/test_distribution.py
index b0fecef2..d2e0c29a 100644
--- a/HARK/tests/test_distribution.py
+++ b/HARK/tests/test_distribution.py
@@ -30,13 +30,32 @@ class DistributionClassTests(unittest.TestCase):
             3.5397367004222002)
 
     def test_Lognormal(self):
+
+        dist = distribution.Lognormal()
+        
+        self.assertEqual(
+            dist.draw(1)[0],
+            5.836039190663969)
+
+        dist.draw(100)
+        dist.reset()
+
         self.assertEqual(
-            distribution.Lognormal().draw(1)[0],
+            dist.draw(1)[0],
             5.836039190663969)
 
     def test_Normal(self):
+        dist = distribution.Normal()
+        
+        self.assertEqual(
+            dist.draw(1)[0],
+            1.764052345967664)
+
+        dist.draw(100)
+        dist.reset()
+
         self.assertEqual(
-            distribution.Normal().draw(1)[0],
+            dist.draw(1)[0],
             1.764052345967664)
 
     def test_Weibull(self):