Merge pull request #355 from ariddell/feature/issue-341-transformed-data-rng

Use random seed in transformed data rng

Author: ariddell

pystan/model.py

@@ -468,8 +468,8 @@ def optimizing(self, data=None, seed=None,
             raise ValueError("Algorithm must be one of {}".format(algorithms))
         if data is None:
             data = {}
-
-        fit = self.fit_class(data)
+        seed = pystan.misc._check_seed(seed)
+        fit = self.fit_class(data, seed)
 
         m_pars = fit._get_param_names()
         p_dims = fit._get_param_dims()
@@ -486,8 +486,6 @@ def optimizing(self, data=None, seed=None,
                 not isinstance(init, string_types):
             raise ValueError("Wrong specification of initial values.")
 
-        seed = pystan.misc._check_seed(seed)
-
         stan_args = dict(init=init,
                          seed=seed,
                          method="optim",
@@ -670,7 +668,8 @@ def sampling(self, data=None, pars=None, chains=4, iter=2000,
         if algorithm not in algorithms:
             raise ValueError("Algorithm must be one of {}".format(algorithms))
 
-        fit = self.fit_class(data)
+        seed = pystan.misc._check_seed(seed)
+        fit = self.fit_class(data, seed)
 
         m_pars = fit._get_param_names()
         p_dims = fit._get_param_dims()
@@ -848,7 +847,8 @@ def vb(self, data=None, pars=None, iter=10000,
         algorithm = "meanfield" if algorithm is None else algorithm
         if algorithm not in algorithms:
             raise ValueError("Algorithm must be one of {}".format(algorithms))
-        fit = self.fit_class(data)
+        seed = pystan.misc._check_seed(seed)
+        fit = self.fit_class(data, seed)
         m_pars = fit._get_param_names()
         p_dims = fit._get_param_dims()
 
@@ -860,8 +860,6 @@ def vb(self, data=None, pars=None, iter=10000,
                 not isinstance(init, string_types):
             raise ValueError("Wrong specification of initial values.")
 
-        seed = pystan.misc._check_seed(seed)
-
         stan_args = dict(iter=iter,
                          init=init,
                          chain_id=1,

pystan/stan_fit.hpp

@@ -1371,9 +1371,9 @@ namespace pystan {
       return true;
     }
 
-    stan_fit(vars_r_t& vars_r, vars_i_t& vars_i) :
+    stan_fit(vars_r_t& vars_r, vars_i_t& vars_i, unsigned int random_seed) :
       data_(vars_r, vars_i),
-      model_(data_, &std::cout),
+      model_(data_, random_seed, &std::cout),
       base_rng(static_cast<boost::uint32_t>(std::time(0))),
       names_(get_param_names(model_)),
       dims_(get_param_dims(model_)),

pystan/stan_fit.pxd

@@ -40,7 +40,7 @@ cdef extern from "stan_fit.hpp" namespace "pystan":
         VARIATIONAL = 4
 
     cdef cppclass stan_fit[M, R]:
-        stan_fit(vars_r_t& vars_r, vars_i_t& vars_i) except +
+        stan_fit(vars_r_t& vars_r, vars_i_t& vars_i, uint random_seed) except +
         bool update_param_oi(vector[string] pars)
         vector[double] unconstrain_pars(vars_r_t& vars_r, vars_i_t& vars_i)
         vector[double] constrain_pars(vector[double]& params_r) except +

pystan/stanfit4model.pyx

@@ -388,7 +388,7 @@ def _call_sampler(data, args, pars_oi=None):
     cdef stan_fit[$model_cppname, ecuyer1988] *fitptr
     cdef vars_r_t vars_r = _dict_to_vars_r(data_r)
     cdef vars_i_t vars_i = _dict_to_vars_i(data_i)
-    fitptr = new stan_fit[$model_cppname, ecuyer1988](vars_r, vars_i)
+    fitptr = new stan_fit[$model_cppname, ecuyer1988](vars_r, vars_i, argsptr.random_seed)
     if not fitptr:
         raise MemoryError("Couldn't allocate space for stan_fit.")
     # Implementation note: there is an extra stan_fit instance associated
@@ -434,6 +434,7 @@ cdef class StanFit4Model:
 
     # attributes populated by methods of StanModel
     cdef public data  # dict or OrderedDict
+    cdef public random_seed
     cdef public dict sim
     cdef public model_name
     cdef public model_pars
@@ -448,18 +449,20 @@ cdef class StanFit4Model:
         # __cinit__ must be callable with no arguments for unpickling
         cdef vars_r_t vars_r
         cdef vars_i_t vars_i
-        if len(args) == 1:
-            data = args[0]
+        if len(args) == 2:
+            data, random_seed = args
             data_r, data_i = pystan.misc._split_data(data)
             # NB: dictionary keys must be byte strings
             vars_r = _dict_to_vars_r(data_r)
             vars_i = _dict_to_vars_i(data_i)
-            self.thisptr = new stan_fit[$model_cppname, ecuyer1988](vars_r, vars_i)
+            # TODO: the random seed needs to be known by StanFit4Model
+            self.thisptr = new stan_fit[$model_cppname, ecuyer1988](vars_r, vars_i, <unsigned int> random_seed)
             if not self.thisptr:
                 raise MemoryError("Couldn't allocate space for stan_fit.")
 
-    def __init__(self, data):
+    def __init__(self, data, random_seed):
         self.data = data
+        self.random_seed = random_seed
 
     def __dealloc__(self):
         del self.thisptr
@@ -480,7 +483,7 @@ cdef class StanFit4Model:
                "The relevant StanModel instance must be pickled along with this fit object.\n"
                "When unpickling the StanModel must be unpickled first.")
         warnings.warn(msg)
-        return (StanFit4Model, (self.data,), self.__getstate__(), None, None)
+        return (StanFit4Model, (self.data, self.random_seed), self.__getstate__(), None, None)
 
     # public methods
 

pystan/tests/test_generated_quantities_seed.py

@@ -0,0 +1,53 @@
+from collections import OrderedDict
+import gc
+import os
+import tempfile
+import unittest
+
+import numpy as np
+
+import pystan
+from pystan._compat import PY2
+
+
+class TestGeneratedQuantitiesSeed(unittest.TestCase):
+    """Verify that the RNG in the transformed data block uses the overall seed.
+
+    See https://github.com/stan-dev/stan/issues/2241
+
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        model_code = """
+            data {
+              int<lower=0> N;
+            }
+            transformed data {
+              vector[N] y;
+              for (n in 1:N)
+                y[n] = normal_rng(0, 1);
+            }
+            parameters {
+              real mu;
+              real<lower = 0> sigma;
+            }
+            model {
+              y ~ normal(mu, sigma);
+            }
+            generated quantities {
+              real mean_y = mean(y);
+              real sd_y = sd(y);
+            }
+        """
+        cls.model = pystan.StanModel(model_code=model_code, verbose=True)
+
+    def test_generated_quantities_seed(self):
+        fit1 = self.model.sampling(data={'N': 1000}, iter=10, seed=123)
+        extr1 = fit1.extract()
+        fit2 = self.model.sampling(data={'N': 1000}, iter=10, seed=123)
+        extr2 = fit2.extract()
+        self.assertTrue((extr1['mean_y'] == extr2['mean_y']).all())
+        fit3 = self.model.sampling(data={'N': 1000}, iter=10, seed=456)
+        extr3 = fit3.extract()
+        self.assertFalse((extr1['mean_y'] == extr3['mean_y']).all())