RNG in transformed data should use seed

[bob-carpenter]

Summary:

The RNG in the transformed data block is insensitive to seed. It should use the overall seed, but not the chain ID.

Thanks to Jim Savage for reporting.

Reproducible Steps:

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);
}

Current Output:

In RStan, the following two calls produce exactly the same fit:

> fit <- stan("foo.stan", data=list(N=1000), seed=123)

> fit <- stan("foo.stan", data=list(N=1000), seed=848475)

as does any seed that's included.

Expected Output:

Different results based on seed.

Additional Information:

The RNG is probably being used before the seed is set in the command interface. We should check how this works post-refactor as I'm testing it in the pre-refactor current-release version.

Current Version:

v2.14.0

[mitzimorris]

this is a problem for all of the interfaces. the generated model contains both 2- and 3-arg constructors, e.g., using the above stand code in file "m2241.stan", the file "m2241.hpp",
has the following two constructors:

    m2241_model(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        typedef boost::ecuyer1988 rng_t;
        rng_t base_rng(0);  // 0 seed default
        ctor_body(context__, base_rng, pstream__);
    }

    template <class RNG>
    m2241_model(stan::io::var_context& context__,
        RNG& base_rng__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, base_rng__, pstream__);
    }

in cmdstan, in file src/cmdstan/command.hpp line 102, the model is instantiated using the 2-arg constructor. the fully instantiated constructor is passed into the calls to stan::services::{diagnose,optimize,sample}

if a random seed is specified, this should be changed to:

boost::ecuyer1988 rng(random_seed);
Model model(data_var_context, rng, &std::cout);

[mitzimorris]

made above change cmdstan.
given model, as above with addition of print(y) statement, i.e.:

transformed data {
  vector[N] y;
  for (n in 1:N)
    y[n] = normal_rng(0, 1);
  print(y);
}

running model with seeds 12345 and 45678 produces different initial values for vector y using modified cmdstan, whereas current version of cmdstan produces identical values, as expected, given that 2-arg constructor always instantiates rng with seed 0.

[mitzimorris]

closing this as stan issue - will add issue for interfaces.

[ariddell]

Could we change the stan::services functions so that the interfaces don't have to keep track of random seeds in two places? It's actually rather nice (and in keeping with the spirit of the refactor) to use the 2-arg constructor. The interfaces shouldn't have to worry about what RNG is being used (beyond the chain id and the integer random seed), right?

[mitzimorris]

would it make sense to make the arg to the 3-arg constructor be the integer random seed instead of a boost rng object?
the interfaces are responsible for instantiating the model object which is passed in to the stan::services functions.

[ariddell]

Now that I look at the code (stan::services::sample code), I worry there might be another concern. Does the random sampling need to vary across chains?

transformed data {
  vector[N] y;
  for (n in 1:N)
    y[n] = normal_rng(0, 1);
}

What's the expected behavior here given a fixed seed and multiple chains?

[ariddell]

@mitzimorris , having the third argument be the random seed sounds like a good idea to me. I had originally hoped one could construct the stan_model instance inside the stan::services function call.

I imagined changing this

      template <class Model>
      int hmc_nuts_diag_e(Model& model, stan::io::var_context& init,
                          unsigned int random_seed, unsigned int chain,
                          ...) {
        boost::ecuyer1988 rng = util::create_rng(random_seed, chain);

into something approximately like this:

      template <class Model>
      int hmc_nuts_diag_e(stan::io::var_context& data, stan::io::var_context& init,
                          unsigned int random_seed, unsigned int chain,
                          ...) {
        Model model = Model(data, base_rng, anything_else);
        boost::ecuyer1988 rng = util::create_rng(random_seed, chain);

but I have a suspicion that this might exceed the limits of what is allowed with templates.

Your solution solves the problem just as well though. It does require the interfaces to remember to pass the random seed in two places -- but that's better than having to worry about the base_rng class AND the random seed in two places.

[mitzimorris]

Does the random sampling need to vary across chains?

we're talking about transformed data block - what do you mean by "sampling"? this is executed once, as part of model's constructor.

if the transformed data should be randomized differently across chains, then the interface will need to instantiate the model on a per-chain basis, passing in both the seed and chain.

having the third argument be the random seed sounds like a good idea to me

ran this by @bob-carpenter who says it's a bad idea - suggestion withdrawn.

[bob-carpenter]

I believe we had decided that the randomization should vary across chains; we could change our minds here.

In some ways, it makes sense not to vary transformed data RNGs by chain ID---no idea how we'd diagnose convergence otherwise. In other ways, the point of it being random is that you want it to vary.

Either way, we want to reuse the RNG for transformed data and then for sampling.

What Allen suggests will work. You can indeed do this:

template <typename M>
...
  stan::io::var_context data;
  boost::ecuyer1988 rng = util::create_rng(random_seed, chain);
  // (1) advance rng by seed here to have varying inits per chain
  M model(data, rng);
  // (2) advance rng by seed here to have same init per chain

The way templates work is that you get real types substituted for all template parameters before compliation, so that will look just like any other type constructor call.

Either way, this preserves the encapsulation of the model not knowing about chain IDs.

[sakrejda]

I thought the idea was that we promise deterministic behavior if you set chain ID and we don't promise anything else. If somebody wants to do a simulation study with repeat data used in multiple places they can just pass it in as data.

[bob-carpenter]

The issue is convergence monitoring. I'm just trying to think like Andrew here. If we have different random transformed data across seeds, then we can't monitor convergence with R-hat and our n_eff will be weird. So I think we should actually have the same RNG for the transformed data, independent of seed. I don't think we're promising anything particularly clearly. Maybe we should.

[sakrejda]

I like your rationale Bob.

…

On Tue, May 16, 2017, 12:39 PM Bob Carpenter ***@***.***> wrote: The issue is convergence monitoring. I'm just trying to think like Andrew here. If we have different random transformed data across seeds, then we can't monitor convergence with R-hat and our n_eff will be weird. So I think we should actually have the same RNG for the transformed data, independent of seed. I don't think we're promising anything particularly clearly. Maybe we should. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2241 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAfA6TqFo-MMaWXtcBjVip_rI9SY-5bvks5r6dFIgaJpZM4MU5zw> .

[bob-carpenter]

I meant independent of chain ID, not independent of seed! So hard to talk about all this precisely with the chain ID acting as a kind of sub-seed (it's used to step the RNG over a bunch of draws).

[mitzimorris]

independent of chain ID, not independent of seed

in which case, why not pass in seed instead of RNG to model constructor?

[betanalpha]

There are two use cases here. First, having constant transformed data to ensure that we can use convergence monitoring over parallel chains. Second, facilitating the sample different data sets and fit to determine ensemble behavior. Both are valid but being able to do both would require keeping track of multiple seeds and being clever about default behavior. I suggest that we use a common seed right now and at least get reproducible behavior and then worry about generalizing to include both behaviors later. A long answer that reduces to: “what Bob said”.

…

On May 16, 2017, at 12:39 PM, Bob Carpenter ***@***.***> wrote: The issue is convergence monitoring. I'm just trying to think like Andrew here. If we have different random transformed data across seeds, then we can't monitor convergence with R-hat and our n_eff will be weird. So I think we should actually have the same RNG for the transformed data, independent of seed. I don't think we're promising anything particularly clearly. Maybe we should. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2241 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABdNlj66840ru-5usIwOQEGG7M2tsp3Mks5r6dFJgaJpZM4MU5zw>.

[bob-carpenter]

@mitzimorris Yes, we just need the ID. We can advance all the other RNGs at least one notch (a gazillion draws) so they won't overlap.

[mitzimorris]

fixed in #2313