From 15ba923e29ea7c5eb2276fc0a526e9eda2f883b2 Mon Sep 17 00:00:00 2001
From: aloctavodia <aloctavodia@gmail.com>
Date: Sat, 7 Jan 2017 12:56:25 -0300
Subject: [PATCH 01/27] initialize chains from estimated posterior samples

---
 pymc3/sampling.py | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index 6c80b49472..ac68d2d680 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -100,10 +100,14 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
         automatically (defaults to None).
     init : str {'advi', 'advi_map', 'map', 'nuts', None}
         Initialization method to use.
-        * advi : Run ADVI to estimate posterior mean and diagonal covariance matrix.
+        * advi : Run ADVI to estimate starting points and diagonal covariance
+        matrix. If njobs > 1 it will sample starting points from the estimated
+        posterior, otherwise it will use the estimated posterior mean.
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map : Use the MAP as starting point.
-        * nuts : Run NUTS and estimate posterior mean and covariance matrix.
+        * nuts : Run NUTS to estimate starting points and covariance matrix. If
+        njobs > 1 it will sample starting points from the estimated posterior,
+        otherwise it will use the estimated posterior mean.
         * None : Do not initialize.
     n_init : int
         Number of iterations of initializer
@@ -146,7 +150,7 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
     if step is None and init is not None and pm.model.all_continuous(model.vars):
         # By default, use NUTS sampler
         pm._log.info('Auto-assigning NUTS sampler...')
-        start_, step = init_nuts(init=init, n_init=n_init, model=model)
+        start_, step = init_nuts(init=init, njobs=njobs, n_init=n_init, model=model)
         if start is None:
             start = start_
     else:
@@ -393,7 +397,7 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see
     return {k: np.asarray(v) for k, v in ppc.items()}
 
 
-def init_nuts(init='advi', n_init=500000, model=None, **kwargs):
+def init_nuts(init='advi', njobs=1, n_init=500000, model=None, **kwargs):
     """Initialize and sample from posterior of a continuous model.
 
     This is a convenience function. NUTS convergence and sampling speed is extremely
@@ -409,6 +413,8 @@ def init_nuts(init='advi', n_init=500000, model=None, **kwargs):
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map : Use the MAP as starting point.
         * nuts : Run NUTS and estimate posterior mean and covariance matrix.
+    njobs : int
+        Number of parallel jobs to start. 
     n_init : int
         Number of iterations of initializer
         If 'advi', number of iterations, if 'metropolis', number of draws.
@@ -432,7 +438,10 @@ def init_nuts(init='advi', n_init=500000, model=None, **kwargs):
 
     if init == 'advi':
         v_params = pm.variational.advi(n=n_init)
-        start = pm.variational.sample_vp(v_params, 1, progressbar=False, hide_transformed=False)[0]
+        if njobs > 1:
+            start = pm.variational.sample_vp(v_params, njobs, progressbar=False, hide_transformed=False)
+        else:
+            start = v_params.means
         cov = np.power(model.dict_to_array(v_params.stds), 2)
     elif init == 'advi_map':
         start = pm.find_MAP()
@@ -441,12 +450,13 @@ def init_nuts(init='advi', n_init=500000, model=None, **kwargs):
     elif init == 'map':
         start = pm.find_MAP()
         cov = pm.find_hessian(point=start)
-
     elif init == 'nuts':
-        init_trace = pm.sample(step=pm.NUTS(), draws=n_init)
-        cov = pm.trace_cov(init_trace[n_init//2:])
-
-        start = {varname: np.mean(init_trace[varname]) for varname in init_trace.varnames}
+        init_trace = pm.sample(step=pm.NUTS(), draws=n_init)[n_init//2:]
+        cov = np.atleast_1d(pm.trace_cov(init_trace))
+        if njobs > 1:
+            start = np.random.choice(init_trace, njobs)
+        else:
+            start = {varname: np.mean(init_trace[varname]) for varname in init_trace.varnames}
     else:
         raise NotImplemented('Initializer {} is not supported.'.format(init))
 

From e13d1868073888bc64d80efd8a68e54f0b5d1c97 Mon Sep 17 00:00:00 2001
From: aloctavodia <aloctavodia@gmail.com>
Date: Sun, 8 Jan 2017 14:42:05 -0300
Subject: [PATCH 02/27] add random_seed to fix starting points

---
 pymc3/sampling.py | 31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index ac68d2d680..20c3a07cda 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -146,11 +146,12 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
     MultiTrace object with access to sampling values
     """
     model = modelcontext(model)
+    
 
     if step is None and init is not None and pm.model.all_continuous(model.vars):
         # By default, use NUTS sampler
         pm._log.info('Auto-assigning NUTS sampler...')
-        start_, step = init_nuts(init=init, njobs=njobs, n_init=n_init, model=model)
+        start_, step = init_nuts(init=init, njobs=njobs, n_init=n_init, model=model, random_seed=random_seed)
         if start is None:
             start = start_
     else:
@@ -397,7 +398,8 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see
     return {k: np.asarray(v) for k, v in ppc.items()}
 
 
-def init_nuts(init='advi', njobs=1, n_init=500000, model=None, **kwargs):
+def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
+            random_seed=-1, **kwargs):
     """Initialize and sample from posterior of a continuous model.
 
     This is a convenience function. NUTS convergence and sampling speed is extremely
@@ -436,27 +438,30 @@ def init_nuts(init='advi', njobs=1, n_init=500000, model=None, **kwargs):
 
     pm._log.info('Initializing NUTS using {}...'.format(init))
 
+    random_seed = int(np.atleast_1d(random_seed)[0])
+
     if init == 'advi':
-        v_params = pm.variational.advi(n=n_init)
-        if njobs > 1:
-            start = pm.variational.sample_vp(v_params, njobs, progressbar=False, hide_transformed=False)
-        else:
-            start = v_params.means
+        v_params = pm.variational.advi(n=n_init, random_seed=random_seed)
+        start = pm.variational.sample_vp(v_params, njobs, progressbar=False,
+        hide_transformed=False, random_seed=random_seed)
+        if njobs == 1:
+            start = start[0]
         cov = np.power(model.dict_to_array(v_params.stds), 2)
     elif init == 'advi_map':
         start = pm.find_MAP()
-        v_params = pm.variational.advi(n=n_init, start=start)
+        v_params = pm.variational.advi(n=n_init, start=start,
+        random_seed=random_seed)
         cov = np.power(model.dict_to_array(v_params.stds), 2)
     elif init == 'map':
         start = pm.find_MAP()
         cov = pm.find_hessian(point=start)
     elif init == 'nuts':
-        init_trace = pm.sample(step=pm.NUTS(), draws=n_init)[n_init//2:]
+        init_trace = pm.sample(step=pm.NUTS(), draws=n_init,
+        random_seed=random_seed)[n_init//2:]
         cov = np.atleast_1d(pm.trace_cov(init_trace))
-        if njobs > 1:
-            start = np.random.choice(init_trace, njobs)
-        else:
-            start = {varname: np.mean(init_trace[varname]) for varname in init_trace.varnames}
+        start = np.random.choice(init_trace, njobs)
+        if njobs == 1:
+            start = start[0]
     else:
         raise NotImplemented('Initializer {} is not supported.'.format(init))
 

From df34cafc2fa878f34278b601c73cc662ae166e3e Mon Sep 17 00:00:00 2001
From: fredcallaway <fredc@llaway.com>
Date: Mon, 9 Jan 2017 14:11:10 -0800
Subject: [PATCH 03/27] allow one-sided patsy formula in linear_component

---
 pymc3/glm/glm.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pymc3/glm/glm.py b/pymc3/glm/glm.py
index 88c121354c..208248f198 100644
--- a/pymc3/glm/glm.py
+++ b/pymc3/glm/glm.py
@@ -42,6 +42,7 @@ class linear_component(namedtuple('Estimate', 'y_est,coeffs')):
     ----------
     formula : str
         Patsy linear model descriptor.
+        E.g. 'y ~ a + b' or simply 'a + b'
     data : array
         Labeled array (e.g. pandas DataFrame, recarray).
     priors : dict
@@ -86,7 +87,11 @@ def __new__(cls, formula, data, priors=None,
             priors = defaultdict(None)
 
         # Build patsy design matrix and get regressor names.
-        _, dmatrix = patsy.dmatrices(formula, data)
+        # Formula can be of form 'c ~ a + b' or 'a + b'.
+        if '~' in formula:
+            _, dmatrix = patsy.dmatrices(formula, data)
+        else:
+            dmatrix = patsy.dmatrix(formula, data)
         reg_names = dmatrix.design_info.column_names
 
         if init_vals is None:

From 078addd788c9dcef0c4a1e8a6087a043a47ddf06 Mon Sep 17 00:00:00 2001
From: fredcallaway <fredc@llaway.com>
Date: Mon, 9 Jan 2017 18:33:46 -0800
Subject: [PATCH 04/27] fix bug inspecting patsy formula

---
 pymc3/glm/glm.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/pymc3/glm/glm.py b/pymc3/glm/glm.py
index 208248f198..ae8d029b74 100644
--- a/pymc3/glm/glm.py
+++ b/pymc3/glm/glm.py
@@ -87,11 +87,10 @@ def __new__(cls, formula, data, priors=None,
             priors = defaultdict(None)
 
         # Build patsy design matrix and get regressor names.
-        # Formula can be of form 'c ~ a + b' or 'a + b'.
-        if '~' in formula:
-            _, dmatrix = patsy.dmatrices(formula, data)
-        else:
-            dmatrix = patsy.dmatrix(formula, data)
+        try:
+            _, dmatrix = patsy.dmatrices(formula, data)  # 'c ~ a + b'
+        except patsy.PatsyError:
+            dmatrix = patsy.dmatrix(formula, data)  # 'a + b'
         reg_names = dmatrix.design_info.column_names
 
         if init_vals is None:

From 0f5eba37271a576f0b776db1a7c10bf68254a61e Mon Sep 17 00:00:00 2001
From: gladomat <glad@posteo.de>
Date: Wed, 18 Jan 2017 11:42:56 +0100
Subject: [PATCH 05/27] Added output for pointwise predictive accuracy in waic.

---
 pymc3/stats.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pymc3/stats.py b/pymc3/stats.py
index 5bea76ea8c..7a95622cc3 100644
--- a/pymc3/stats.py
+++ b/pymc3/stats.py
@@ -105,7 +105,7 @@ def log_post_trace(trace, model):
     return np.vstack([obs.logp_elemwise(pt) for obs in model.observed_RVs] for pt in trace)
 
 
-def waic(trace, model=None, n_eff=False):
+def waic(trace, model=None, n_eff=False, pointwise=False):
     """
     Calculate the widely available information criterion, its standard error
     and the effective number of parameters of the samples in trace from model.
@@ -121,6 +121,9 @@ def waic(trace, model=None, n_eff=False):
     n_eff: bool
         if True the effective number parameters will be returned.
         Default False
+    pointwise: bool
+        if True the pointwise predictive accuracy will be returned.
+        Default False
 
 
     Returns
@@ -152,6 +155,8 @@ def waic(trace, model=None, n_eff=False):
 
     if n_eff:
         return waic, waic_se, p_waic
+    elif pointwise:
+        return waic, waic_se, waic_i, p_waic
     else:
         return waic, waic_se
 

From a2baf670e8859f017efce8dd4e845ecc398fcd23 Mon Sep 17 00:00:00 2001
From: aloctavodia <aloctavodia@gmail.com>
Date: Wed, 18 Jan 2017 10:40:48 -0300
Subject: [PATCH 06/27] fix indentation

---
 pymc3/sampling.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index 20c3a07cda..69931ca21c 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -399,7 +399,7 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see
 
 
 def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
-            random_seed=-1, **kwargs):
+              random_seed=-1, **kwargs):
     """Initialize and sample from posterior of a continuous model.
 
     This is a convenience function. NUTS convergence and sampling speed is extremely
@@ -443,21 +443,22 @@ def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
     if init == 'advi':
         v_params = pm.variational.advi(n=n_init, random_seed=random_seed)
         start = pm.variational.sample_vp(v_params, njobs, progressbar=False,
-        hide_transformed=False, random_seed=random_seed)
+                                         hide_transformed=False, 
+                                         random_seed=random_seed)
         if njobs == 1:
             start = start[0]
         cov = np.power(model.dict_to_array(v_params.stds), 2)
     elif init == 'advi_map':
         start = pm.find_MAP()
-        v_params = pm.variational.advi(n=n_init, start=start,
-        random_seed=random_seed)
+        v_params = pm.variational.advi(n=n_init, start=start, 
+                                       random_seed=random_seed)
         cov = np.power(model.dict_to_array(v_params.stds), 2)
     elif init == 'map':
         start = pm.find_MAP()
         cov = pm.find_hessian(point=start)
     elif init == 'nuts':
         init_trace = pm.sample(step=pm.NUTS(), draws=n_init,
-        random_seed=random_seed)[n_init//2:]
+                               random_seed=random_seed)[n_init // 2:]
         cov = np.atleast_1d(pm.trace_cov(init_trace))
         start = np.random.choice(init_trace, njobs)
         if njobs == 1:

From da67ffa00cae00afd5bdce0f90892c56de35e4d3 Mon Sep 17 00:00:00 2001
From: Chris Fonnesbeck <fonnesbeck@gmail.com>
Date: Wed, 18 Jan 2017 10:20:56 -0600
Subject: [PATCH 07/27] Changed case of init arguments in sample to upper case,
 but kept backward-compatibility with lower case

---
 pymc3/examples/arbitrary_stochastic.py |  2 +-
 pymc3/examples/factor_potential.py     |  2 +-
 pymc3/sampling.py                      | 31 +++++++++++++++-----------
 pymc3/tests/test_examples.py           |  6 ++---
 pymc3/tests/test_glm.py                |  2 +-
 pymc3/tests/test_models_linear.py      |  8 +++----
 pymc3/tests/test_plots.py              |  4 ++--
 7 files changed, 30 insertions(+), 25 deletions(-)

diff --git a/pymc3/examples/arbitrary_stochastic.py b/pymc3/examples/arbitrary_stochastic.py
index 9739c56dac..a342a453e5 100644
--- a/pymc3/examples/arbitrary_stochastic.py
+++ b/pymc3/examples/arbitrary_stochastic.py
@@ -20,7 +20,7 @@ def run(n_samples=3000):
     start = model.test_point
     h = pm.find_hessian(start, model=model)
     step = pm.Metropolis(model.vars, h, blocked=True, model=model)
-    trace = pm.sample(n_samples, step, start, model=model)
+    trace = pm.sample(n_samples, step=step, start=start, model=model)
     return trace
 
 if __name__ == "__main__":
diff --git a/pymc3/examples/factor_potential.py b/pymc3/examples/factor_potential.py
index 17b1d6c743..a189d729e2 100644
--- a/pymc3/examples/factor_potential.py
+++ b/pymc3/examples/factor_potential.py
@@ -13,7 +13,7 @@ def run(n=3000):
     if n == "short":
         n = 50
     with model:
-        trace = sample(n, step, start)
+        trace = sample(n, step=step, start=start)
 
 if __name__ == '__main__':
     run()
diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index 69931ca21c..95f5660e89 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -98,20 +98,20 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
         A step function or collection of functions. If no step methods are
         specified, or are partially specified, they will be assigned
         automatically (defaults to None).
-    init : str {'advi', 'advi_map', 'map', 'nuts', None}
+    init : str {'ADVI', 'ADVI_MAP', 'MAP', 'NUTS', None}
         Initialization method to use.
-        * advi : Run ADVI to estimate starting points and diagonal covariance
+        * ADVI : Run ADVI to estimate starting points and diagonal covariance
         matrix. If njobs > 1 it will sample starting points from the estimated
         posterior, otherwise it will use the estimated posterior mean.
-        * advi_map: Initialize ADVI with MAP and use MAP as starting point.
-        * map : Use the MAP as starting point.
-        * nuts : Run NUTS to estimate starting points and covariance matrix. If
+        * ADVI_MAP: Initialize ADVI with MAP and use MAP as starting point.
+        * MAP : Use the MAP as starting point.
+        * NUTS : Run NUTS to estimate starting points and covariance matrix. If
         njobs > 1 it will sample starting points from the estimated posterior,
         otherwise it will use the estimated posterior mean.
         * None : Do not initialize.
     n_init : int
         Number of iterations of initializer
-        If 'advi', number of iterations, if 'nuts', number of draws.
+        If 'ADVI', number of iterations, if 'nuts', number of draws.
     start : dict
         Starting point in parameter space (or partial point)
         Defaults to trace.point(-1)) if there is a trace provided and
@@ -147,6 +147,8 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None,
     """
     model = modelcontext(model)
     
+    if init is not None:
+        init = init.lower()
 
     if step is None and init is not None and pm.model.all_continuous(model.vars):
         # By default, use NUTS sampler
@@ -398,7 +400,7 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see
     return {k: np.asarray(v) for k, v in ppc.items()}
 
 
-def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
+def init_nuts(init='ADVI', njobs=1, n_init=500000, model=None,
               random_seed=-1, **kwargs):
     """Initialize and sample from posterior of a continuous model.
 
@@ -409,17 +411,17 @@ def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
 
     Parameters
     ----------
-    init : str {'advi', 'advi_map', 'map', 'nuts'}
+    init : str {'ADVI', 'ADVI_MAP', 'MAP', 'NUTS'}
         Initialization method to use.
-        * advi : Run ADVI to estimate posterior mean and diagonal covariance matrix.
-        * advi_map: Initialize ADVI with MAP and use MAP as starting point.
-        * map : Use the MAP as starting point.
-        * nuts : Run NUTS and estimate posterior mean and covariance matrix.
+        * ADVI : Run ADVI to estimate posterior mean and diagonal covariance matrix.
+        * ADVI_MAP: Initialize ADVI with MAP and use MAP as starting point.
+        * MAP : Use the MAP as starting point.
+        * NUTS : Run NUTS and estimate posterior mean and covariance matrix.
     njobs : int
         Number of parallel jobs to start. 
     n_init : int
         Number of iterations of initializer
-        If 'advi', number of iterations, if 'metropolis', number of draws.
+        If 'ADVI', number of iterations, if 'metropolis', number of draws.
     model : Model (optional if in `with` context)
     **kwargs : keyword arguments
         Extra keyword arguments are forwarded to pymc3.NUTS.
@@ -439,6 +441,9 @@ def init_nuts(init='advi', njobs=1, n_init=500000, model=None,
     pm._log.info('Initializing NUTS using {}...'.format(init))
 
     random_seed = int(np.atleast_1d(random_seed)[0])
+    
+    if init is not None:
+        init = init.lower()
 
     if init == 'advi':
         v_params = pm.variational.advi(n=n_init, random_seed=random_seed)
diff --git a/pymc3/tests/test_examples.py b/pymc3/tests/test_examples.py
index 2ce8f1a793..c2da6c1c2b 100644
--- a/pymc3/tests/test_examples.py
+++ b/pymc3/tests/test_examples.py
@@ -54,7 +54,7 @@ def test_run(self):
             h = H(start)
 
             step = pm.HamiltonianMC(model.vars, h)
-            pm.sample(50, step, start)
+            pm.sample(50, step=step, start=start)
 
 
 class TestARM12_6(SeededTest):
@@ -88,7 +88,7 @@ def too_slow(self):
             start = pm.find_MAP(start=start,
                                 vars=[model['groupmean'], model['sd_interval_'], model['floor_m']])
             step = pm.NUTS(model.vars, scaling=start)
-            pm.sample(50, step, start)
+            pm.sample(50, step=step, start=start)
 
 
 class TestARM12_6Uranium(SeededTest):
@@ -129,7 +129,7 @@ def too_slow(self):
             h = np.diag(H(start))
 
             step = pm.HamiltonianMC(model.vars, h)
-            pm.sample(50, step, start)
+            pm.sample(50, step=step, start=start)
 
 
 def build_disaster_model(masked=False):
diff --git a/pymc3/tests/test_glm.py b/pymc3/tests/test_glm.py
index 259ffdb48b..b4cc107ed7 100644
--- a/pymc3/tests/test_glm.py
+++ b/pymc3/tests/test_glm.py
@@ -34,7 +34,7 @@ def test_linear_component(self):
             Normal('y_obs', mu=y_est, sd=sigma, observed=self.y_linear)
             start = find_MAP(vars=[sigma])
             step = Slice(model.vars)
-            trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed)
+            trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed)
 
             self.assertAlmostEqual(np.mean(trace['Intercept']), self.intercept, 1)
             self.assertAlmostEqual(np.mean(trace['x']), self.slope, 1)
diff --git a/pymc3/tests/test_models_linear.py b/pymc3/tests/test_models_linear.py
index fac66fc0de..b0a00953e7 100644
--- a/pymc3/tests/test_models_linear.py
+++ b/pymc3/tests/test_models_linear.py
@@ -44,7 +44,7 @@ def test_linear_component(self):
             Normal('y_obs', mu=lm.y_est, sd=sigma, observed=self.y_linear)  # yields y_obs
             start = find_MAP(vars=[sigma])
             step = Slice(model.vars)
-            trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed)
+            trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed)
 
             self.assertAlmostEqual(np.mean(trace['lm_Intercept']), self.intercept, 1)
             self.assertAlmostEqual(np.mean(trace['lm_x0']), self.slope, 1)
@@ -58,7 +58,7 @@ def test_linear_component_from_formula(self):
             Normal('y_obs', mu=lm.y_est, sd=sigma, observed=self.y_linear)
             start = find_MAP(vars=[sigma])
             step = Slice(model.vars)
-            trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed)
+            trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed)
 
             self.assertAlmostEqual(np.mean(trace['Intercept']), self.intercept, 1)
             self.assertAlmostEqual(np.mean(trace['x']), self.slope, 1)
@@ -79,7 +79,7 @@ def test_glm(self):
             )
             start = find_MAP()
             step = Slice(model.vars)
-            trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed)
+            trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed)
             self.assertAlmostEqual(np.mean(trace['glm_Intercept']), self.intercept, 1)
             self.assertAlmostEqual(np.mean(trace['glm_x0']), self.slope, 1)
             self.assertAlmostEqual(np.mean(trace['glm_sd']), self.sd, 1)
@@ -91,7 +91,7 @@ def test_glm_from_formula(self):
             Glm.from_formula('y ~ x', self.data_linear, name=NAME)
             start = find_MAP()
             step = Slice(model.vars)
-            trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed)
+            trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed)
 
             self.assertAlmostEqual(np.mean(trace['%s_Intercept' % NAME]), self.intercept, 1)
             self.assertAlmostEqual(np.mean(trace['%s_x' % NAME]), self.slope, 1)
diff --git a/pymc3/tests/test_plots.py b/pymc3/tests/test_plots.py
index 05dc8d5a46..dc2ea22e44 100644
--- a/pymc3/tests/test_plots.py
+++ b/pymc3/tests/test_plots.py
@@ -20,7 +20,7 @@ def test_plots():
         start = model.test_point
         h = find_hessian(start)
         step = Metropolis(model.vars, h)
-        trace = sample(3000, step, start)
+        trace = sample(3000, step=step, start=start)
 
         traceplot(trace)
         forestplot(trace)
@@ -34,7 +34,7 @@ def test_plots_multidimensional():
     with model as model:
         h = np.diag(find_hessian(start))
         step = Metropolis(model.vars, h)
-        trace = sample(3000, step, start)
+        trace = sample(3000, step=step, start=start)
 
         traceplot(trace)
         plot_posterior(trace)

From 016e80784d4a4e51291f0f6e1b160224e3e169c1 Mon Sep 17 00:00:00 2001
From: Chris Fonnesbeck <fonnesbeck@gmail.com>
Date: Wed, 18 Jan 2017 12:29:00 -0600
Subject: [PATCH 08/27] Changed positional arguments to keyword arguments in
 tests

---
 pymc3/tests/test_diagnostics.py | 14 +++++++-------
 pymc3/tests/test_plots.py       |  2 +-
 pymc3/tests/test_sampling.py    |  2 +-
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/pymc3/tests/test_diagnostics.py b/pymc3/tests/test_diagnostics.py
index d5e6987385..770ca2ce7d 100644
--- a/pymc3/tests/test_diagnostics.py
+++ b/pymc3/tests/test_diagnostics.py
@@ -20,9 +20,9 @@ def get_ptrace(self, n_samples):
             # Run sampler
             step1 = Slice([model.early_mean_log_, model.late_mean_log_])
             step2 = Metropolis([model.switchpoint])
-            start = {'early_mean': 2., 'late_mean': 3., 'switchpoint': 90}
-            ptrace = sample(n_samples, [step1, step2], start, njobs=2, progressbar=False,
-                            random_seed=[1, 4])
+            start = {'early_mean': 7., 'late_mean': 1., 'switchpoint': 90}
+            ptrace = sample(n_samples, step=[step1, step2], start=start, njobs=2, 
+                    progressbar=False, random_seed=[1, 4])
         return ptrace
 
     def test_good(self):
@@ -53,7 +53,7 @@ def test_right_shape_python_float(self, shape=None, test_shape=None):
             # start sampling at the MAP
             start = find_MAP()
             step = NUTS(scaling=start)
-            ptrace = sample(n_samples, step, start,
+            ptrace = sample(n_samples, step=step, start=start,
                             njobs=n_jobs, random_seed=42)
 
         rhat = gelman_rubin(ptrace)['x']
@@ -92,7 +92,7 @@ def get_switchpoint(self, n_samples):
             # Run sampler
             step1 = Slice([model.early_mean_log_, model.late_mean_log_])
             step2 = Metropolis([model.switchpoint])
-            trace = sample(n_samples, [step1, step2], progressbar=False, random_seed=1)
+            trace = sample(n_samples, step=[step1, step2], progressbar=False, random_seed=1)
         return trace['switchpoint']
 
     def test_geweke_negative(self):
@@ -153,7 +153,7 @@ def test_effective_n(self):
             # start sampling at the MAP
             start = find_MAP()
             step = NUTS(scaling=start)
-            ptrace = sample(n_samples, step, start,
+            ptrace = sample(n_samples, step=step, start=start,
                             njobs=n_jobs, random_seed=42)
 
         n_effective = effective_n(ptrace)['x']
@@ -174,7 +174,7 @@ def test_effective_n_right_shape_python_float(self,
             # start sampling at the MAP
             start = find_MAP()
             step = NUTS(scaling=start)
-            ptrace = sample(n_samples, step, start,
+            ptrace = sample(n_samples, step=step, start=start,
                             njobs=n_jobs, random_seed=42)
 
         n_effective = effective_n(ptrace)['x']
diff --git a/pymc3/tests/test_plots.py b/pymc3/tests/test_plots.py
index dc2ea22e44..a0d23464da 100644
--- a/pymc3/tests/test_plots.py
+++ b/pymc3/tests/test_plots.py
@@ -47,7 +47,7 @@ def test_multichain_plots():
         step1 = Slice([model.early_mean_log_, model.late_mean_log_])
         step2 = Metropolis([model.switchpoint])
         start = {'early_mean': 2., 'late_mean': 3., 'switchpoint': 50}
-        ptrace = sample(1000, [step1, step2], start, njobs=2)
+        ptrace = sample(1000, step=[step1, step2], start=start, njobs=2)
 
     forestplot(ptrace, varnames=['early_mean', 'late_mean'])
     autocorrplot(ptrace, varnames=['switchpoint'])
diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py
index 5c81be566c..7d85d2284d 100644
--- a/pymc3/tests/test_sampling.py
+++ b/pymc3/tests/test_sampling.py
@@ -59,7 +59,7 @@ def test_sample(self):
         with self.model:
             for njobs in test_njobs:
                 for steps in [1, 10, 300]:
-                    pm.sample(steps, self.step, {}, None, njobs=njobs, random_seed=self.random_seed)
+                    pm.sample(steps, step=self.step, njobs=njobs, random_seed=self.random_seed)
 
     def test_sample_init(self):
         with self.model:

From faab2e1f6b2679909d15bec76c2a6cd803657ac9 Mon Sep 17 00:00:00 2001
From: Chris Fonnesbeck <fonnesbeck@gmail.com>
Date: Wed, 18 Jan 2017 16:29:35 -0600
Subject: [PATCH 09/27] Fix failures in test_examples and test_diagnostics

---
 pymc3/tests/test_diagnostics.py | 4 ++--
 pymc3/tests/test_examples.py    | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pymc3/tests/test_diagnostics.py b/pymc3/tests/test_diagnostics.py
index 770ca2ce7d..d85738012d 100644
--- a/pymc3/tests/test_diagnostics.py
+++ b/pymc3/tests/test_diagnostics.py
@@ -20,9 +20,9 @@ def get_ptrace(self, n_samples):
             # Run sampler
             step1 = Slice([model.early_mean_log_, model.late_mean_log_])
             step2 = Metropolis([model.switchpoint])
-            start = {'early_mean': 7., 'late_mean': 1., 'switchpoint': 90}
+            start = {'early_mean': 7., 'late_mean': 1., 'switchpoint': 100}
             ptrace = sample(n_samples, step=[step1, step2], start=start, njobs=2, 
-                    progressbar=False, random_seed=[1, 4])
+                    progressbar=False, random_seed=[20090425, 19700903])
         return ptrace
 
     def test_good(self):
diff --git a/pymc3/tests/test_examples.py b/pymc3/tests/test_examples.py
index c2da6c1c2b..d071826cce 100644
--- a/pymc3/tests/test_examples.py
+++ b/pymc3/tests/test_examples.py
@@ -263,7 +263,7 @@ def test_run(self):
             start = {'psi': 0.5, 'z': (self.y > 0).astype(int), 'theta': 5}
             step_one = pm.Metropolis([model.theta_interval_, model.psi_logodds_])
             step_two = pm.BinaryMetropolis([model.z])
-            pm.sample(50, [step_one, step_two], start)
+            pm.sample(50, step=[step_one, step_two], start=start)
 
 
 class TestRSV(SeededTest):

From 20aa34690805f1906a6e9a524d1382bc5d035fb2 Mon Sep 17 00:00:00 2001
From: Chris Fonnesbeck <fonnesbeck@gmail.com>
Date: Thu, 17 Nov 2016 13:43:51 -0500
Subject: [PATCH 10/27] Added mode argument to several step methods and advi to
 allow mode setting (e.g. nanguardmode) for Theano functions

---
 pymc3/step_methods/metropolis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pymc3/step_methods/metropolis.py b/pymc3/step_methods/metropolis.py
index c2ab7260ca..4fcbdc7866 100644
--- a/pymc3/step_methods/metropolis.py
+++ b/pymc3/step_methods/metropolis.py
@@ -433,6 +433,6 @@ def delta_logp(logp, vars, shared):
 
     logp1 = pm.CallableTensor(logp0)(inarray1)
 
-    f = theano.function([inarray1, inarray0], logp1 - logp0)
+    f = theano.function([inarray1, inarray0], logp1 - logp0, mode=self.mode)
     f.trust_input = True
     return f

From 12c556a13b913dbaf1478ad41f1012ae9fcae84f Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 19:21:20 +0300
Subject: [PATCH 11/27] Created Generator Op with simple Test

---
 pymc3/data.py               | 34 +++++++++++++++++++++++++++++++++-
 pymc3/tests/test_theanof.py | 18 ++++++++++++++++++
 pymc3/theanof.py            | 32 ++++++++++++++++++++++++++++++++
 3 files changed, 83 insertions(+), 1 deletion(-)
 create mode 100644 pymc3/tests/test_theanof.py

diff --git a/pymc3/data.py b/pymc3/data.py
index 944ce55458..bd0a6e5735 100644
--- a/pymc3/data.py
+++ b/pymc3/data.py
@@ -1,7 +1,10 @@
+import itertools
 import pkgutil
 import io
 
-__all__ = ['get_data_file']
+import theano.tensor as tt
+
+__all__ = ['get_data_file', 'DataGenerator']
 
 
 def get_data_file(pkg, path):
@@ -19,3 +22,32 @@ def get_data_file(pkg, path):
     """
 
     return io.BytesIO(pkgutil.get_data(pkg, path))
+
+
+class DataGenerator(object):
+    """
+    Helper class that helps to infer data type of generator with looking
+    at the first item, preserving the order of the resulting generator
+    """
+    def __init__(self, iterable):
+        if hasattr(iterable, '__len__'):
+            generator = itertools.cycle(iterable)
+        else:
+            generator = iter(iterable)
+        self.test_value = next(generator)
+        self.gen = itertools.chain([self.test_value], generator)
+        self.tensortype = tt.TensorType(
+            self.test_value.dtype,
+            ((False, ) * self.test_value.ndim))
+
+    def __next__(self):
+        return next(self.gen)
+
+    def __iter__(self):
+        return self
+
+    def __eq__(self, other):
+        return id(self) == id(other)
+
+    def __hash__(self):
+        return hash(id(self))
diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
new file mode 100644
index 0000000000..c7fa50ab3e
--- /dev/null
+++ b/pymc3/tests/test_theanof.py
@@ -0,0 +1,18 @@
+import unittest
+import numpy as np
+import theano
+from ..theanof import DataGenerator, GeneratorOp
+
+
+class TestGenerator(unittest.TestCase):
+    def test_basic(self):
+        def integers():
+            i = 0
+            while True:
+                yield np.float32(i)
+                i += 1
+        generator = DataGenerator(integers())
+        gop = GeneratorOp(generator)()
+        f = theano.function([], gop)
+        self.assertEqual(f(), np.float32(0))
+        self.assertEqual(f(), np.float32(1))
diff --git a/pymc3/theanof.py b/pymc3/theanof.py
index 56bad3fb0e..8e8d4c6af8 100644
--- a/pymc3/theanof.py
+++ b/pymc3/theanof.py
@@ -3,8 +3,12 @@
 from .vartypes import typefilter, continuous_types
 from theano import theano, scalar, tensor as tt
 from theano.gof.graph import inputs
+from theano.tensor import TensorVariable
+from theano.gof import Op, Apply
+from theano.configparser import change_flags
 from .memoize import memoize
 from .blocking import ArrayOrdering
+from .data import DataGenerator
 
 __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars',
            'cont_inputs', 'floatX', 'jacobian',
@@ -242,3 +246,31 @@ def __call__(self, input):
 
 scalar_identity = IdentityOp(scalar.upgrade_to_float, name='scalar_identity')
 identity = tt.Elemwise(scalar_identity, name='identity')
+
+
+class GeneratorOp(Op):
+    __props__ = ('generator',)
+
+    def __init__(self, generator):
+        if not isinstance(generator, DataGenerator):
+            raise ValueError('pymc3 requires special DataGenerator for consistency, '
+                             'wrap your generator with it')
+        super(GeneratorOp, self).__init__()
+        self.generator = generator
+        self.itypes = []
+        self.otypes = [self.generator.tensortype]
+
+    def perform(self, node, inputs, output_storage, params=None):
+        output_storage[0][0] = self.generator.__next__()
+
+    def do_constant_folding(self, node):
+        return False
+
+    @change_flags(compute_test_value='off')
+    def __call__(self, *args, **kwargs):
+        rval = super(GeneratorOp, self).__call__(*args, **kwargs)
+        rval.tag.test_value = self.generator.test_value
+        return rval
+
+
+

From b842b3353f17dd630b8f8803cc5c436292e8c8d5 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 19:37:25 +0300
Subject: [PATCH 12/27] added ndim test

---
 pymc3/tests/test_theanof.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
index c7fa50ab3e..31a66a943c 100644
--- a/pymc3/tests/test_theanof.py
+++ b/pymc3/tests/test_theanof.py
@@ -1,3 +1,4 @@
+import itertools
 import unittest
 import numpy as np
 import theano
@@ -16,3 +17,18 @@ def integers():
         f = theano.function([], gop)
         self.assertEqual(f(), np.float32(0))
         self.assertEqual(f(), np.float32(1))
+
+    def test_ndim(self):
+        for ndim in range(10):
+            def integers():
+                i = 0
+                while True:
+                    yield np.ones((2,) * ndim) * i
+                    i += 1
+            res = list(itertools.islice(integers(), 0, 2))
+            generator = DataGenerator(integers())
+            gop = GeneratorOp(generator)()
+            f = theano.function([], gop)
+            self.assertEqual(ndim, res[0].ndim)
+            np.testing.assert_equal(f(), res[0])
+            np.testing.assert_equal(f(), res[1])

From 52b4f3a34096a19b4a02b138c696ba2ecb56acc6 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 19:39:13 +0300
Subject: [PATCH 13/27] updated test

---
 pymc3/tests/test_theanof.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
index 31a66a943c..758b5a3173 100644
--- a/pymc3/tests/test_theanof.py
+++ b/pymc3/tests/test_theanof.py
@@ -17,6 +17,9 @@ def integers():
         f = theano.function([], gop)
         self.assertEqual(f(), np.float32(0))
         self.assertEqual(f(), np.float32(1))
+        for i in range(2, 100):
+            f()
+        self.assertEqual(f(), np.float32(100))
 
     def test_ndim(self):
         for ndim in range(10):

From dced530b93673d406bbd04967060e4c1d69c9acf Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 19:43:11 +0300
Subject: [PATCH 14/27] updated test, added test value check

---
 pymc3/tests/test_theanof.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
index 758b5a3173..86c32416a7 100644
--- a/pymc3/tests/test_theanof.py
+++ b/pymc3/tests/test_theanof.py
@@ -14,6 +14,7 @@ def integers():
                 i += 1
         generator = DataGenerator(integers())
         gop = GeneratorOp(generator)()
+        self.assertEqual(gop.tag.test_value, np.float32(0))
         f = theano.function([], gop)
         self.assertEqual(f(), np.float32(0))
         self.assertEqual(f(), np.float32(1))

From 76940e79292c10ba27e4fe57367c37cdaeea000c Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 19:49:26 +0300
Subject: [PATCH 15/27] added test for replacing generator with shared variable

---
 pymc3/tests/test_theanof.py | 37 +++++++++++++++++++++++++------------
 1 file changed, 25 insertions(+), 12 deletions(-)

diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
index 86c32416a7..7907b98da6 100644
--- a/pymc3/tests/test_theanof.py
+++ b/pymc3/tests/test_theanof.py
@@ -5,13 +5,22 @@
 from ..theanof import DataGenerator, GeneratorOp
 
 
+def integers():
+    i = 0
+    while True:
+        yield np.float32(i)
+        i += 1
+
+
+def integers_ndim(ndim):
+    i = 0
+    while True:
+        yield np.ones((2,) * ndim) * i
+        i += 1
+
+
 class TestGenerator(unittest.TestCase):
     def test_basic(self):
-        def integers():
-            i = 0
-            while True:
-                yield np.float32(i)
-                i += 1
         generator = DataGenerator(integers())
         gop = GeneratorOp(generator)()
         self.assertEqual(gop.tag.test_value, np.float32(0))
@@ -24,15 +33,19 @@ def integers():
 
     def test_ndim(self):
         for ndim in range(10):
-            def integers():
-                i = 0
-                while True:
-                    yield np.ones((2,) * ndim) * i
-                    i += 1
-            res = list(itertools.islice(integers(), 0, 2))
-            generator = DataGenerator(integers())
+            res = list(itertools.islice(integers_ndim(ndim), 0, 2))
+            generator = DataGenerator(integers_ndim(ndim))
             gop = GeneratorOp(generator)()
             f = theano.function([], gop)
             self.assertEqual(ndim, res[0].ndim)
             np.testing.assert_equal(f(), res[0])
             np.testing.assert_equal(f(), res[1])
+
+    def test_cloning_available(self):
+        generator = DataGenerator(integers())
+        gop = GeneratorOp(generator)()
+        res = gop ** 2
+        shared = theano.shared(np.float32(10))
+        res1 = theano.clone(res, {gop: shared})
+        f = theano.function([], res1)
+        self.assertEqual(f(), np.float32(100))

From a4a42075032c89f51f7eb40601b3a90a6f7ec717 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 20:01:17 +0300
Subject: [PATCH 16/27] added shortcut for generator op

---
 pymc3/theanof.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/pymc3/theanof.py b/pymc3/theanof.py
index 8e8d4c6af8..666d66220c 100644
--- a/pymc3/theanof.py
+++ b/pymc3/theanof.py
@@ -4,7 +4,7 @@
 from theano import theano, scalar, tensor as tt
 from theano.gof.graph import inputs
 from theano.tensor import TensorVariable
-from theano.gof import Op, Apply
+from theano.gof import Op
 from theano.configparser import change_flags
 from .memoize import memoize
 from .blocking import ArrayOrdering
@@ -13,7 +13,7 @@
 __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars',
            'cont_inputs', 'floatX', 'jacobian',
            'CallableTensor', 'join_nonshared_inputs',
-           'make_shared_replacements']
+           'make_shared_replacements', 'generator']
 
 
 def inputvars(a):
@@ -251,17 +251,19 @@ def __call__(self, input):
 class GeneratorOp(Op):
     __props__ = ('generator',)
 
-    def __init__(self, generator):
-        if not isinstance(generator, DataGenerator):
-            raise ValueError('pymc3 requires special DataGenerator for consistency, '
-                             'wrap your generator with it')
+    def __init__(self, gen):
+        if not isinstance(gen, DataGenerator):
+            gen = DataGenerator(gen)
         super(GeneratorOp, self).__init__()
-        self.generator = generator
+        self.generator = gen
         self.itypes = []
         self.otypes = [self.generator.tensortype]
 
     def perform(self, node, inputs, output_storage, params=None):
-        output_storage[0][0] = self.generator.__next__()
+        try:
+            output_storage[0][0] = self.generator.__next__()
+        except StopIteration:
+            output_storage[0][0] = np.nan
 
     def do_constant_folding(self, node):
         return False
@@ -273,4 +275,5 @@ def __call__(self, *args, **kwargs):
         return rval
 
 
-
+def generator(gen):
+    return GeneratorOp(gen)()

From 13d5fb3f88594608baeed6ce6bb16f5de9b7425b Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 20:07:49 +0300
Subject: [PATCH 17/27] refactored test

---
 pymc3/tests/test_theanof.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py
index 7907b98da6..ebc7308781 100644
--- a/pymc3/tests/test_theanof.py
+++ b/pymc3/tests/test_theanof.py
@@ -2,7 +2,7 @@
 import unittest
 import numpy as np
 import theano
-from ..theanof import DataGenerator, GeneratorOp
+from ..theanof import DataGenerator, GeneratorOp, generator
 
 
 def integers():
@@ -42,8 +42,7 @@ def test_ndim(self):
             np.testing.assert_equal(f(), res[1])
 
     def test_cloning_available(self):
-        generator = DataGenerator(integers())
-        gop = GeneratorOp(generator)()
+        gop = generator(integers())
         res = gop ** 2
         shared = theano.shared(np.float32(10))
         res1 = theano.clone(res, {gop: shared})

From 7104be56b550999424dbabe7fcab7a624156a2b4 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 20:27:55 +0300
Subject: [PATCH 18/27] added population kwarg (no tests yet)

---
 pymc3/data.py                       |  6 +-----
 pymc3/distributions/distribution.py |  3 ++-
 pymc3/model.py                      | 29 +++++++++++++++++++++--------
 3 files changed, 24 insertions(+), 14 deletions(-)

diff --git a/pymc3/data.py b/pymc3/data.py
index bd0a6e5735..fcac737719 100644
--- a/pymc3/data.py
+++ b/pymc3/data.py
@@ -29,11 +29,7 @@ class DataGenerator(object):
     Helper class that helps to infer data type of generator with looking
     at the first item, preserving the order of the resulting generator
     """
-    def __init__(self, iterable):
-        if hasattr(iterable, '__len__'):
-            generator = itertools.cycle(iterable)
-        else:
-            generator = iter(iterable)
+    def __init__(self, generator):
         self.test_value = next(generator)
         self.gen = itertools.chain([self.test_value], generator)
         self.tensortype = tt.TensorType(
diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index 2ad2243f74..d6ab06d683 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -30,8 +30,9 @@ def __new__(cls, name, *args, **kwargs):
 
         if isinstance(name, string_types):
             data = kwargs.pop('observed', None)
+            population = kwargs.pop('population', None)
             dist = cls.dist(*args, **kwargs)
-            return model.Var(name, dist, data)
+            return model.Var(name, dist, data, population)
         else:
             raise TypeError("Name needs to be a string but got: %s" % name)
 
diff --git a/pymc3/model.py b/pymc3/model.py
index 6017484d06..7c34f6047e 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -8,7 +8,7 @@
 
 import pymc3 as pm
 from .memoize import memoize
-from .theanof import gradient, hessian, inputvars
+from .theanof import gradient, hessian, inputvars, generator
 from .vartypes import typefilter, discrete_types, continuous_types
 from .blocking import DictToArrayBijection, ArrayOrdering
 
@@ -458,7 +458,7 @@ def cont_vars(self):
         """All the continuous variables in the model"""
         return list(typefilter(self.vars, continuous_types))
 
-    def Var(self, name, dist, data=None):
+    def Var(self, name, dist, data=None, population=None):
         """Create and add (un)observed random variable to the model with an
         appropriate prior distribution.
 
@@ -491,7 +491,7 @@ def Var(self, name, dist, data=None):
                 return var
         elif isinstance(data, dict):
             var = MultiObservedRV(name=name, data=data, distribution=dist,
-                                  model=self)
+                                  model=self, population=population)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs += var.missing_values
@@ -500,7 +500,8 @@ def Var(self, name, dist, data=None):
                     self.named_vars[v.name] = v
         else:
             var = ObservedRV(name=name, data=data,
-                             distribution=dist, model=self)
+                             distribution=dist, model=self,
+                             population=population)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs.append(var.missing_values)
@@ -780,6 +781,8 @@ def as_tensor(data, name, model, distribution):
             constant[data.mask.nonzero()], missing_values)
         dataTensor.missing_values = missing_values
         return dataTensor
+    elif hasattr(data, '__next__'):
+        return generator(data)
     else:
         data = tt.as_tensor_variable(data, name=name)
         data.missing_values = None
@@ -792,7 +795,7 @@ class ObservedRV(Factor, TensorVariable):
     """
 
     def __init__(self, type=None, owner=None, index=None, name=None, data=None,
-                 distribution=None, model=None):
+                 distribution=None, model=None, population=None):
         """
         Parameters
         ----------
@@ -814,7 +817,12 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
             data = as_tensor(data, name, model, distribution)
             self.missing_values = data.missing_values
 
-            self.logp_elemwiset = distribution.logp(data)
+            logp_elemwiset = distribution.logp(data)
+            if population is None:
+                coef = tt.as_tensor(1)
+            else:
+                coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+            self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
             self.distribution = distribution
 
@@ -835,7 +843,7 @@ class MultiObservedRV(Factor):
     Potentially partially observed.
     """
 
-    def __init__(self, name, data, distribution, model):
+    def __init__(self, name, data, distribution, model, population=None):
         """
         Parameters
         ----------
@@ -851,7 +859,12 @@ def __init__(self, name, data, distribution, model):
 
         self.missing_values = [datum.missing_values for datum in self.data.values()
                                if datum.missing_values is not None]
-        self.logp_elemwiset = distribution.logp(**self.data)
+        logp_elemwiset = distribution.logp(**self.data)
+        if population is None:
+            coef = tt.as_tensor(1)
+        else:
+            coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+        self.logp_elemwiset = logp_elemwiset * coef
         self.model = model
         self.distribution = distribution
 

From a530e7d2970f0e91b53cd70115651296b7a84830 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Wed, 11 Jan 2017 21:20:23 +0300
Subject: [PATCH 19/27] added population kwarg for free var(autoencoder case)

---
 pymc3/model.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/pymc3/model.py b/pymc3/model.py
index 7c34f6047e..c437c55e75 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -477,11 +477,13 @@ def Var(self, name, dist, data=None, population=None):
         name = self.name_for(name)
         if data is None:
             if getattr(dist, "transform", None) is None:
-                var = FreeRV(name=name, distribution=dist, model=self)
+                var = FreeRV(name=name, distribution=dist, model=self,
+                             population=population)
                 self.free_RVs.append(var)
             else:
                 var = TransformedRV(name=name, distribution=dist, model=self,
-                                    transform=dist.transform)
+                                    transform=dist.transform,
+                                    population=population)
                 pm._log.debug('Applied {transform}-transform to {name}'
                               ' and added transformed {orig_name} to model.'.format(
                                 transform=dist.transform.name,
@@ -718,7 +720,7 @@ class FreeRV(Factor, TensorVariable):
     """Unobserved random variable that a model is specified in terms of."""
 
     def __init__(self, type=None, owner=None, index=None, name=None,
-                 distribution=None, model=None):
+                 distribution=None, model=None, population=None):
         """
         Parameters
         ----------
@@ -737,7 +739,12 @@ def __init__(self, type=None, owner=None, index=None, name=None,
             self.distribution = distribution
             self.tag.test_value = np.ones(
                 distribution.shape, distribution.dtype) * distribution.default()
-            self.logp_elemwiset = distribution.logp(self)
+            logp_elemwiset = distribution.logp(self)
+            if population is None:
+                coef = tt.as_tensor(1)
+            else:
+                coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+            self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
 
             incorporate_methods(source=distribution, destination=self,
@@ -909,7 +916,8 @@ def Potential(name, var, model=None):
 class TransformedRV(TensorVariable):
 
     def __init__(self, type=None, owner=None, index=None, name=None,
-                 distribution=None, model=None, transform=None):
+                 distribution=None, model=None, transform=None,
+                 population=None):
         """
         Parameters
         ----------
@@ -929,7 +937,7 @@ def __init__(self, type=None, owner=None, index=None, name=None,
 
             transformed_name = "{}_{}_".format(name, transform.name)
             self.transformed = model.Var(
-                transformed_name, transform.apply(distribution))
+                transformed_name, transform.apply(distribution), population=population)
 
             normalRV = transform.backward(self.transformed)
 

From babd88c7927b38ab5c6018b817706becffc341eb Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Thu, 12 Jan 2017 13:16:28 +0300
Subject: [PATCH 20/27] Revert "Added mode argument to several step methods and
 advi to allow mode setting (e.g. nanguardmode) for Theano functions"

This reverts commit 6bdb53a069c6d569fadcbc36a93489c0ddb334a3.
---
 pymc3/step_methods/metropolis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pymc3/step_methods/metropolis.py b/pymc3/step_methods/metropolis.py
index 4fcbdc7866..c2ab7260ca 100644
--- a/pymc3/step_methods/metropolis.py
+++ b/pymc3/step_methods/metropolis.py
@@ -433,6 +433,6 @@ def delta_logp(logp, vars, shared):
 
     logp1 = pm.CallableTensor(logp0)(inarray1)
 
-    f = theano.function([inarray1, inarray0], logp1 - logp0, mode=self.mode)
+    f = theano.function([inarray1, inarray0], logp1 - logp0)
     f.trust_input = True
     return f

From 50d2abfeda3009f36ed2ba62b32d91e8cdb4f36a Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Sat, 14 Jan 2017 13:25:48 +0300
Subject: [PATCH 21/27] add docstring to generator Op

---
 pymc3/theanof.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/pymc3/theanof.py b/pymc3/theanof.py
index 666d66220c..eeddacf5ce 100644
--- a/pymc3/theanof.py
+++ b/pymc3/theanof.py
@@ -249,6 +249,17 @@ def __call__(self, input):
 
 
 class GeneratorOp(Op):
+    """
+    Generaror Op is designed for storing python generators
+    inside theano graph. The main limitation is generator itself.
+
+    There are some important cases when generator becomes exhausted
+        - not endless generator is passed
+        - exception is raised while `generator.__next__` is performed
+            Note: it is dangerous in simple python generators, but ok in
+            custom class based generators with explicit state
+        - data type on each iteration should be the same
+    """
     __props__ = ('generator',)
 
     def __init__(self, gen):
@@ -276,4 +287,5 @@ def __call__(self, *args, **kwargs):
 
 
 def generator(gen):
+    """shortcut for `GeneratorOp`"""
     return GeneratorOp(gen)()

From 2d6f23874464d20ba379abdc6140e15b018968ad Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Sat, 14 Jan 2017 13:27:54 +0300
Subject: [PATCH 22/27] rename population -> total_size

---
 pymc3/distributions/distribution.py |  4 ++--
 pymc3/model.py                      | 32 ++++++++++++++---------------
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index d6ab06d683..cae9740740 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -30,9 +30,9 @@ def __new__(cls, name, *args, **kwargs):
 
         if isinstance(name, string_types):
             data = kwargs.pop('observed', None)
-            population = kwargs.pop('population', None)
+            total_size = kwargs.pop('total_size', None)
             dist = cls.dist(*args, **kwargs)
-            return model.Var(name, dist, data, population)
+            return model.Var(name, dist, data, total_size)
         else:
             raise TypeError("Name needs to be a string but got: %s" % name)
 
diff --git a/pymc3/model.py b/pymc3/model.py
index c437c55e75..46ce029b61 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -458,7 +458,7 @@ def cont_vars(self):
         """All the continuous variables in the model"""
         return list(typefilter(self.vars, continuous_types))
 
-    def Var(self, name, dist, data=None, population=None):
+    def Var(self, name, dist, data=None, total_size=None):
         """Create and add (un)observed random variable to the model with an
         appropriate prior distribution.
 
@@ -478,12 +478,12 @@ def Var(self, name, dist, data=None, population=None):
         if data is None:
             if getattr(dist, "transform", None) is None:
                 var = FreeRV(name=name, distribution=dist, model=self,
-                             population=population)
+                             total_size=total_size)
                 self.free_RVs.append(var)
             else:
                 var = TransformedRV(name=name, distribution=dist, model=self,
                                     transform=dist.transform,
-                                    population=population)
+                                    total_size=total_size)
                 pm._log.debug('Applied {transform}-transform to {name}'
                               ' and added transformed {orig_name} to model.'.format(
                                 transform=dist.transform.name,
@@ -493,7 +493,7 @@ def Var(self, name, dist, data=None, population=None):
                 return var
         elif isinstance(data, dict):
             var = MultiObservedRV(name=name, data=data, distribution=dist,
-                                  model=self, population=population)
+                                  model=self, total_size=total_size)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs += var.missing_values
@@ -503,7 +503,7 @@ def Var(self, name, dist, data=None, population=None):
         else:
             var = ObservedRV(name=name, data=data,
                              distribution=dist, model=self,
-                             population=population)
+                             total_size=total_size)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs.append(var.missing_values)
@@ -720,7 +720,7 @@ class FreeRV(Factor, TensorVariable):
     """Unobserved random variable that a model is specified in terms of."""
 
     def __init__(self, type=None, owner=None, index=None, name=None,
-                 distribution=None, model=None, population=None):
+                 distribution=None, model=None, total_size=None):
         """
         Parameters
         ----------
@@ -740,10 +740,10 @@ def __init__(self, type=None, owner=None, index=None, name=None,
             self.tag.test_value = np.ones(
                 distribution.shape, distribution.dtype) * distribution.default()
             logp_elemwiset = distribution.logp(self)
-            if population is None:
+            if total_size is None:
                 coef = tt.as_tensor(1)
             else:
-                coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+                coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
             self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
 
@@ -802,7 +802,7 @@ class ObservedRV(Factor, TensorVariable):
     """
 
     def __init__(self, type=None, owner=None, index=None, name=None, data=None,
-                 distribution=None, model=None, population=None):
+                 distribution=None, model=None, total_size=None):
         """
         Parameters
         ----------
@@ -825,10 +825,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
             self.missing_values = data.missing_values
 
             logp_elemwiset = distribution.logp(data)
-            if population is None:
+            if total_size is None:
                 coef = tt.as_tensor(1)
             else:
-                coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+                coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
             self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
             self.distribution = distribution
@@ -850,7 +850,7 @@ class MultiObservedRV(Factor):
     Potentially partially observed.
     """
 
-    def __init__(self, name, data, distribution, model, population=None):
+    def __init__(self, name, data, distribution, model, total_size=None):
         """
         Parameters
         ----------
@@ -867,10 +867,10 @@ def __init__(self, name, data, distribution, model, population=None):
         self.missing_values = [datum.missing_values for datum in self.data.values()
                                if datum.missing_values is not None]
         logp_elemwiset = distribution.logp(**self.data)
-        if population is None:
+        if total_size is None:
             coef = tt.as_tensor(1)
         else:
-            coef = tt.as_tensor(population) / logp_elemwiset.shape[0]
+            coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
         self.logp_elemwiset = logp_elemwiset * coef
         self.model = model
         self.distribution = distribution
@@ -917,7 +917,7 @@ class TransformedRV(TensorVariable):
 
     def __init__(self, type=None, owner=None, index=None, name=None,
                  distribution=None, model=None, transform=None,
-                 population=None):
+                 total_size=None):
         """
         Parameters
         ----------
@@ -937,7 +937,7 @@ def __init__(self, type=None, owner=None, index=None, name=None,
 
             transformed_name = "{}_{}_".format(name, transform.name)
             self.transformed = model.Var(
-                transformed_name, transform.apply(distribution), population=population)
+                transformed_name, transform.apply(distribution), total_size=total_size)
 
             normalRV = transform.backward(self.transformed)
 

From 9fdf17e91ca22955669521b25f8d43be4b959a88 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Sat, 14 Jan 2017 13:40:53 +0300
Subject: [PATCH 23/27] update docstrings in model

---
 pymc3/model.py | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/pymc3/model.py b/pymc3/model.py
index 46ce029b61..bdfaaf0c15 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -469,6 +469,8 @@ def Var(self, name, dist, data=None, total_size=None):
         data : array_like (optional)
            If data is provided, the variable is observed. If None,
            the variable is unobserved.
+        total_size : scalar
+            upscales logp of variable with :math:`coef = total_size/var.shape[0]`
 
         Returns
         -------
@@ -728,7 +730,10 @@ def __init__(self, type=None, owner=None, index=None, name=None,
         owner : theano owner (optional)
         name : str
         distribution : Distribution
-        model : Model"""
+        model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
+        """
         if type is None:
             type = distribution.type
         super(FreeRV, self).__init__(type, owner, index, name)
@@ -811,6 +816,8 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
         name : str
         distribution : Distribution
         model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
         """
         from .distributions import TensorType
         if type is None:
@@ -859,6 +866,8 @@ def __init__(self, name, data, distribution, model, total_size=None):
         name : str
         distribution : Distribution
         model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
         """
         self.name = name
         self.data = {name: as_tensor(data, name, model, distribution)
@@ -924,10 +933,12 @@ def __init__(self, type=None, owner=None, index=None, name=None,
 
         type : theano type (optional)
         owner : theano owner (optional)
-
         name : str
         distribution : Distribution
-        model : Model"""
+        model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
+        """
         if type is None:
             type = distribution.type
         super(TransformedRV, self).__init__(type, owner, index, name)

From 0ba59bee3e0fdae265837e710904c62943dace9b Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Sat, 14 Jan 2017 13:41:18 +0300
Subject: [PATCH 24/27] fix typo in `as_tensor` function

---
 pymc3/model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pymc3/model.py b/pymc3/model.py
index bdfaaf0c15..56241971fe 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -772,6 +772,8 @@ def pandas_to_array(data):
         return data
     elif isinstance(data, theano.gof.graph.Variable):
         return data
+    elif hasattr(data, '__next__'):
+        return generator(data)
     else:
         return np.asarray(data)
 
@@ -793,8 +795,6 @@ def as_tensor(data, name, model, distribution):
             constant[data.mask.nonzero()], missing_values)
         dataTensor.missing_values = missing_values
         return dataTensor
-    elif hasattr(data, '__next__'):
-        return generator(data)
     else:
         data = tt.as_tensor_variable(data, name=name)
         data.missing_values = None

From 1630a7602bac59c24d749e3db513d29e1dc83054 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Fri, 20 Jan 2017 01:03:38 +0300
Subject: [PATCH 25/27] add test for density rescaling

---
 pymc3/model.py            | 12 +++++++++---
 pymc3/tests/test_model.py | 12 +++++++++++-
 2 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/pymc3/model.py b/pymc3/model.py
index 56241971fe..46d598b7b3 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -746,8 +746,10 @@ def __init__(self, type=None, owner=None, index=None, name=None,
                 distribution.shape, distribution.dtype) * distribution.default()
             logp_elemwiset = distribution.logp(self)
             if total_size is None:
-                coef = tt.as_tensor(1)
+                coef = tt.constant(1)
             else:
+                assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                                  'needs to be at least 1 dimensional')
                 coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
             self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
@@ -833,8 +835,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
 
             logp_elemwiset = distribution.logp(data)
             if total_size is None:
-                coef = tt.as_tensor(1)
+                coef = tt.constant(1)
             else:
+                assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                                  'needs to be at least 1 dimensional')
                 coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
             self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
@@ -877,8 +881,10 @@ def __init__(self, name, data, distribution, model, total_size=None):
                                if datum.missing_values is not None]
         logp_elemwiset = distribution.logp(**self.data)
         if total_size is None:
-            coef = tt.as_tensor(1)
+            coef = tt.constant(1)
         else:
+            assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                              'needs to be at least 1 dimensional')
             coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
         self.logp_elemwiset = logp_elemwiset * coef
         self.model = model
diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py
index 4eb96e6789..dd86c293a9 100644
--- a/pymc3/tests/test_model.py
+++ b/pymc3/tests/test_model.py
@@ -1,5 +1,5 @@
 import unittest
-import theano.tensor as tt
+from theano import theano, tensor as tt
 import pymc3 as pm
 from pymc3.distributions import HalfCauchy, Normal
 from pymc3 import Potential, Deterministic
@@ -118,3 +118,13 @@ def test_model_root(self):
             self.assertTrue(model is model.root)
             with pm.Model() as sub:
                 self.assertTrue(model is sub.root)
+
+    def test_density_scaling(self):
+        with pm.Model() as model1:
+            Normal('n', observed=[[1]], total_size=1)
+            p1 = theano.function([], model1.logpt)
+
+        with pm.Model() as model2:
+            Normal('n', observed=[[1]], total_size=2)
+            p2 = theano.function([], model2.logpt)
+            self.assertEqual(p1() * 2, p2())

From 91011b5870b9d27886f330118739a1e543a6e552 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Fri, 20 Jan 2017 01:15:33 +0300
Subject: [PATCH 26/27] add test for density with generators

---
 pymc3/tests/test_model.py | 31 ++++++++++++++++++++++++++++++-
 1 file changed, 30 insertions(+), 1 deletion(-)

diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py
index dd86c293a9..42d56e667e 100644
--- a/pymc3/tests/test_model.py
+++ b/pymc3/tests/test_model.py
@@ -1,5 +1,6 @@
 import unittest
 from theano import theano, tensor as tt
+import numpy as np
 import pymc3 as pm
 from pymc3.distributions import HalfCauchy, Normal
 from pymc3 import Potential, Deterministic
@@ -127,4 +128,32 @@ def test_density_scaling(self):
         with pm.Model() as model2:
             Normal('n', observed=[[1]], total_size=2)
             p2 = theano.function([], model2.logpt)
-            self.assertEqual(p1() * 2, p2())
+        self.assertEqual(p1() * 2, p2())
+
+    def test_density_scaling_with_genarator(self):
+        # We have different size generators
+        def gen1():
+            i = 0
+            while True:
+                yield np.ones((10, 100)) * i
+                i += 1
+
+        def gen2():
+            i = 0
+            while True:
+                yield np.ones((20, 100)) * i
+                i += 1
+
+        # We have same size models
+        with pm.Model() as model1:
+            Normal('n', observed=gen1(), total_size=100)
+            p1 = theano.function([], model1.logpt)
+
+        with pm.Model() as model2:
+            Normal('n', observed=gen2(), total_size=100)
+            p2 = theano.function([], model2.logpt)
+
+        # We want densities to be equal
+        for _ in range(10):
+            np.testing.assert_almost_equal(p1(), p2())
+        # Done

From 0f36f06efa3aa5e78d61571ba1bdd4ad1d97fc55 Mon Sep 17 00:00:00 2001
From: Maxim Kochurov <maxim.v.kochurov@gmail.com>
Date: Fri, 20 Jan 2017 13:56:11 +0300
Subject: [PATCH 27/27] fix python2 error with generators

---
 pymc3/data.py    | 2 ++
 pymc3/model.py   | 4 +++-
 pymc3/theanof.py | 2 +-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/pymc3/data.py b/pymc3/data.py
index fcac737719..d40d58ca0a 100644
--- a/pymc3/data.py
+++ b/pymc3/data.py
@@ -39,6 +39,8 @@ def __init__(self, generator):
     def __next__(self):
         return next(self.gen)
 
+    next = __next__
+
     def __iter__(self):
         return self
 
diff --git a/pymc3/model.py b/pymc3/model.py
index 46d598b7b3..0a0d5a04b4 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -1,5 +1,6 @@
 import threading
 import six
+import types
 
 import numpy as np
 import theano
@@ -774,7 +775,8 @@ def pandas_to_array(data):
         return data
     elif isinstance(data, theano.gof.graph.Variable):
         return data
-    elif hasattr(data, '__next__'):
+    elif (hasattr(data, '__next__') or
+          isinstance(data, types.GeneratorType)):
         return generator(data)
     else:
         return np.asarray(data)
diff --git a/pymc3/theanof.py b/pymc3/theanof.py
index eeddacf5ce..6870119f7e 100644
--- a/pymc3/theanof.py
+++ b/pymc3/theanof.py
@@ -272,7 +272,7 @@ def __init__(self, gen):
 
     def perform(self, node, inputs, output_storage, params=None):
         try:
-            output_storage[0][0] = self.generator.__next__()
+            output_storage[0][0] = next(self.generator)
         except StopIteration:
             output_storage[0][0] = np.nan