Possible issue with new RNG? about stan HOT 16 CLOSED

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.

from stan.

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

from stan.

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

from stan.

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.

from stan.

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

from stan.

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

from stan.

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

from stan.

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

from stan.

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?

from stan.

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

from stan.

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

from stan.

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.

from stan.

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!

from stan.

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?)

from stan.

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.

from stan.

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

from stan.