Possible issue with new RNG?

[andrjohns]

Description:

Not sure if this is a bug or just me doing things wrong, but it looks like the new RNG doesn't interact with the Math library in the same way - the return is always constant:

#include <iostream>
#include <stan/services/util/create_rng.hpp>
#include <stan/math.hpp>

int main() {
  using stan::math::normal_rng;
  stan::rng_t rng(0);
  std::cout << "stan::rng_t: "
            << normal_rng(5, 10, rng) << ", "
            << normal_rng(5, 10, rng) << ", "
            << normal_rng(5, 10, rng) << ", "
            << normal_rng(5, 10, rng)
            << std::endl;

  boost::ecuyer1988 rng2(0);
  std::cout << "boost::ecuyer1988: "
            << normal_rng(5, 10, rng2) << ", "
            << normal_rng(5, 10, rng2) << ", "
            << normal_rng(5, 10, rng2) << ", "
            << normal_rng(5, 10, rng2)
            << std::endl;
  return 0;
}

Gives:

stan::rng_t: 5, 5, 5, 5
boost::ecuyer1988: -4.52988, 8.12959, 2.12997, 5.78816

Current Version:

v2.35.0

[andrjohns]

Oh hold on - it looks like the issue is the 0 seed value. If I change the seed to 1 then all is good:

stan::rng_t: 0.0297492, 10.3881, -9.09565, -1.47649
boost::ecuyer1988: -4.52988, 8.12959, 2.12997, 5.78816

[bob-carpenter]

That sounds really dangerous---@mitzimorris and @SteveBronder should know if this is going to be an issue through Stan or CmdStan---we need to check somewhere that seeds aren't zero.

[andrjohns]

RStan will also need to update their handling for exposing RNG functions, since they default to a zero seed (FYI @bgoodri) - we were doing the same in cmdstanr, which is how I found this

[bgoodri]

I think that a seed of zero was working fine in the past.

[andrjohns]

I think that a seed of zero was working fine in the past.

Yep, but it will be a problem with the new RNG in Stan (boost::mixmax instead of boost::ecuyer1988) - see the example at the top of the issue

[SteveBronder]

Should we change it to check for zero as the seed and if it is just +1?

[bob-carpenter]

No, because then a user supplying 0 and a user supplying 1 get the same answer. It'd be better to add 1 to all the inputs, which just means whatever the max is will have to be one less.

[WardBrian]

I can't reproduce up at the cmdstan level - supplying seed 0 and id 0, I still get valid sampling behavior, no repeat draws

@andrjohns do you observe the bad behavior up at the cmdstan level? If so this is possibly an architecture issue or boost bug?

[WardBrian]

No, because then a user supplying 0 and a user supplying 1 get the same answer

Note that this was actually already the case with the previous rng in Stan

[andrjohns]

I can't reproduce up at the cmdstan level - supplying seed 0 and id 0, I still get valid sampling behavior, no repeat draws

@andrjohns do you observe the bad behavior up at the cmdstan level? If so this is possibly an architecture issue or boost bug?

Looks like it's definitely a Boost bug - reproduced on godbolt. I'll file an issue with Boost now

EDIT: upstream issue here: boostorg/random#104

[WardBrian]

The papers on the mixmax generator say the seed needs at least 1 nonzero bit. We have a bunch of extra bits to work with in the end, so we can definitely ensure that

[WardBrian]

Separate from the discussion on what's the best way to turn the user-supplied seed into an initialization for the pRNG, do we have an explanation for why this bug seems to not matter at the integration-test level?

E.g., why can I still get good samples from my model if I specify seed=0 id=0 at the moment.

[andrjohns]

Separate from the discussion on what's the best way to turn the user-supplied seed into an initialization for the pRNG, do we have an explanation for why this bug seems to not matter at the integration-test level?

E.g., why can I still get good samples from my model if I specify seed=0 id=0 at the moment.

You can see it with the unconstrained initial values. If you update the random_var_context to print the inits after generation:

      for (size_t n = 0; n < num_unconstrained_; ++n) {
        unconstrained_params_[n] = unif(rng);
        std::cout << "unconstrained_params_[" << n << "] = "
                  << unconstrained_params_[n] << std::endl;
      }

And update the bernoulli example to have multiple parameters (vector<lower=0, upper=1>[N] theta;)

You can see in the output:

andrew@Andrews-MacBook-Air bernoulli % ./bernoulli sample id=0 data file=bernoulli.data.json random seed=0 output diagnostic_file=bernoulli_diags.json
method = sample (Default)
  sample
    num_samples = 1000 (Default)
    num_warmup = 1000 (Default)
    save_warmup = false (Default)
    thin = 1 (Default)
    adapt
      engaged = true (Default)
      gamma = 0.05 (Default)
      delta = 0.8 (Default)
      kappa = 0.75 (Default)
      t0 = 10 (Default)
      init_buffer = 75 (Default)
      term_buffer = 50 (Default)
      window = 25 (Default)
      save_metric = false (Default)
    algorithm = hmc (Default)
      hmc
        engine = nuts (Default)
          nuts
            max_depth = 10 (Default)
        metric = diag_e (Default)
        metric_file =  (Default)
        stepsize = 1 (Default)
        stepsize_jitter = 0 (Default)
    num_chains = 1 (Default)
id = 0
data
  file = bernoulli.data.json
init = 2 (Default)
random
  seed = 0
output
  file = output.csv (Default)
  diagnostic_file = bernoulli_diags.json
  refresh = 100 (Default)
  sig_figs = -1 (Default)
  profile_file = profile.csv (Default)
  save_cmdstan_config = false (Default)
num_threads = 1 (Default)

unconstrained_params_[0] = -2
unconstrained_params_[1] = -2
unconstrained_params_[2] = -2
unconstrained_params_[3] = -2
unconstrained_params_[4] = -2
unconstrained_params_[5] = -2
unconstrained_params_[6] = -2
unconstrained_params_[7] = -2
unconstrained_params_[8] = -2
unconstrained_params_[9] = -2

Gradient evaluation took 2.6e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
Adjust your expectations accordingly!

[WardBrian]

Let me rephrase: at what point does the randomness "start working" such that I end up with what look like valid samples from the posterior with the bad seed? Presumably if the momentum refresh in HMC is also getting turned into a deterministic value, we wouldn't be getting correct sampling behavior (unless the Bernoulli example is so simple it "works" anyway?)

[nhuurre]

Even if momentum resampling is deterministic, the MCMC trajectory is going to be quite complicated and (I imagine) might well look pseudorandom. Maybe you need more effective samples, or the problems are only visible in higher moments (variance, skewness, etc)
If the underlying rng is stuck then _rng in generated quantities won't work at all.

[WardBrian]

I was running a version of Bernoulli which was also doing a posterior predictive check to generate new Ys, which also looked fine, but presumably that’s because the parameters of the rng were different each iteration