Closed cescalara closed 1 year ago
Very odd. I can confirm I can reproduce the behavior you see.
Investigating now
Hm, inserting an extra rng call in the tranformed data block changes the behavior:
data {
real Nex;
}
transformed data {
real ignored = uniform_rng(0,1);
int N_td;
N_td = poisson_rng(Nex);
}
generated quantities {
int N_out_td;
int N_gq;
N_out_td = N_td;
N_gq = poisson_rng(Nex);
}
I suspect the bug you've found is not in cmdstanpy, but somewhere deeper in Stan. I will transfer this issue when I figure out where it should go!
Looking at the code now, it seems we
1) use a different RNG (with the same seed) inside the model constructor as everywhere else
2) do not supply the chain ID to advance the RNG in the constructor, so all chains would have the same RNG in the transformed data block (we just pass chain id 0
).
Both of those seem like potential issues, but neither really explain whatever is going on here. (Edit: 2 is intentional: https://discourse.mc-stan.org/t/transformed-data-rng-varies-by-chain/569)
I've also confirmed that setting the RNG to any other chain (by editing poison.hpp) fixes the problem, and setting id=0
breaks the RNG for generated quantities as well.
Very odd!
Oh also, it has nothing to do with poisson
- here's the same observation using uniform_rng(0,1)
:
Thanks for looking into this! Happy to help out if I can. I thought it was particularly weird that setting lambda>10
seemed to fix things...
It seems like this issue may be even below Stan, but it's possible we can fix it ourselves.
Here's the equivalent of what we're doing in pure C++ with Boost: https://godbolt.org/z/h98Gbdf9j
If you change out seeds, you'll see that the first value printed is always very high, but the following ones look normal.
Very odd! We're still looking into this but we will definitely make sure something is changed for the next version of Stan! For now, inserting a dummy real = uniform_rng(0,1)
call seems to fix subsequent things.
I've transferred this to the Stan repository because this is most likely where it will be fixed, by editing https://github.com/stan-dev/stan/blob/develop/src/stan/services/util/create_rng.hpp
If you change out seeds, you'll see that the first value printed is always very high, but the following ones look normal.
In Stan or in C++?
Just in C++. The godbolt link above shows what I'm describing using only the boost parts we use in Stan.
Just in C++. The godbolt link above shows what I'm describing using only the boost parts we use in Stan.
Are you sure? I did not see anything peculiar?
If you vary the seed, you'll see things like:
Seed: 233
Program stdout
12
2
2
The first number printed is usually >10 regardless of seed, which is quite unlikely for poission(3). Furthermore, if I change the code so that there is always an rng.discard(1)
call, this behavior vanishes
Hmm, the behavior is even weirder, and it seems to also depend on small seeds.
The original python code used a random number between 1 and 10,000 as a seed. If you make this much bigger (I used seed = np.random.randint(10000, 10000000)
) then the graphs start looking correct.
Ok, now I see it
Here's another godbolt which shows the issue more precisely: https://godbolt.org/z/e9veEn5T5
Could you also file an issue with Boost? I know we'll need to work around this in the short term, but it's probably something they want to know, too, as giving things small seeds by hand is a very common practice. My guess is that it's a bad interaction between the linear congruential pRNG we use and the implementation of the Poisson RNG.
It doesn't seem to be unique to Poisson, here is the same code but with a uniform(0,1) rng: https://godbolt.org/z/zfT579dox The first draw is always close to 1, the second draw appears to be from the correct distribution.
That said, it also isn't universal - beta
and normal
both look correct in both situations.
That makes sense. For our pRNG, there's an underlying 256-bit representation of the state of the RNG. I don't know how state is initialized from a 32-bit uint seed, but my guess is that a small number will somehow fail to get a good 256-bit initialization for the way the pRNGs are implemented for Poisson, uniform, etc.
Also, I had no idea C++11 had a lot of this implemented. For example,
https://en.cppreference.com/w/cpp/numeric/random/poisson_distribution
Boost issue opened.
The easiest work-around for us seems to be just always discarding at least 1 draw when we create the RNG. If we want to maximize backwards compatibility of seeds, we can add 1 to the discard amount if and only if the chain ID is 0, otherwise we can just unconditionally add 1 and state in the release notes that random behavior will have changed this version for seeded runs.
@mitzimorris said this issue may have been introduced relatively recently when we changed the default ID of the first chain from 0 to 1 (to avoid having to advance when not required, but it looks like it was required!).
I think changing the default ID to 1 actually made this bug less common. Now it is only apparently in two scenarios
Sorry, I got that backwards. It changed from 1 to 0.
Also, thanks much for reporting, @cescalara. We're all still marveling at how you managed to discover this bug given how subtle and seed-dependent it is.
Looking through the history, it used to discard DISCARD_STRIDE * (chain - 1)
, until this PR in 2017: https://github.com/stan-dev/stan/pull/2313, when it just discarded chain
as it does now. The multi-chain PR https://github.com/stan-dev/stan/pull/3033 edited this file while in progress but ultimately kept this behavior.
I believe the model class has used chain id 0 since the era of that 2017 PR, so the behavior would have been noticeable as long ago as that.
Until we get a boost bugfix, my vote would be to discard std::max(1, DISCARD_STRIDE * chain)
and announce it in release notes and everywhere else. Once/if boost does a bugfix, we then switch back to what we currently have.
The only drawback (at least that I see) is backward compatibility for anyone using chain ID = 0. But anyone using that would have been using a biased rng and there is no reason to leave a bug in for backward compatiblity.
Noone should be using ecuyer1988 RNG as it is obsolete. The replacement boost::ecuyer1988 -> boost::random::mixmax should solve your issues. MIXMAX is a modern, industrial strength RNG. It has been made the default in major physics Monte-Carlo packages, GEANT4 and CMS-SW software. https://github.com/Geant4/geant4/blob/master/source/externals/clhep/src/MixMaxRng.cc https://github.com/cms-sw/cmssw/pull/21339/commits/e79b4675d9d9907c053e4e1a99089271ad6bdeb4
Along with sequences indistinguishable statistically from true random numbers, attention has been paid to seeding, such that the sequence obtained from seed(1) starts out perfectly randomized and is completely independent from seed(2), seed(3) etc, including the case if you want to test the quality by interleaving the two or more streams.
It appears MIXMAX defines discard as
void discard(boost::uint64_t nsteps) { for(boost::uint64_t j = 0; j < nsteps; ++j) (*this)(); } ///< discard n steps, required in boost::random
We regularly call discard with an argument on the order of 2^50, so this wouldn’t exactly be a drop in replacement for our use case and would require reconsidering how we use the random number generators
@WardBrian , @bob-carpenter
Most RNG work on a single trajectory, where seeding can place the state of the RNG at a "random" point.
The RNG called ecuyer1988 has the period of about 10^18, so advancing on it for 2^50 steps would exhaust the period and create the danger that random numbers simply repeat between different draws. Moreover, it is certain that with a generator like ecuyer1988 will produce identical stream of random numbers from pairs of completely different seeds (with a lag of e.g. 1).
In MIXMAX, the trajectory is astronomically long and the state of seed(k+1) is the state of seed(k) advanced by exactly 2^512 steps. Perhaps if I understood what task you accomplish with discard, I could find some solution for you.
The RNG called ecuyer1988 has the period of about 10^18
Table 32.6 in the Boost documentation claims ecuyer1988
has a period of approximately 2^61
Your comment https://github.com/boostorg/random/issues/92#issuecomment-1446210886 is correct about why we discard, it is for parallel MCMC chains. We currently create each chain with the same seed but discard (2^50 * chain_id) elements.
The behavior you're describing for MIXMAX sounds like something we could reasonably use. If I'm understanding correctly we could then seed each chain with (user_supplied_seed + chain_id)
Summary:
When using
poisson_rng
inside the transformed parameter block, the results do not follow the expected Poisson distribution. The issue only seems to appear for rate values oflambda < 10
, e.g. poisson_rng(9.9). The issue only appears inside the transformed data block, and usingpoisson_rng
inside the generated quantities block is fine.Description:
If I compile and run the Stan code shown below as a simulation (i.e.
fixed_param=True
) via cmdstanpy I get very different output depending on which block the function is called. When callingpoisson_rng
in the transformed data block, the output doesn't make sense. I also add the python code I use and a plot demonstrating this.Stan code:
poisson.stan
python code:
Figure:![stan_poisson_db](https://user-images.githubusercontent.com/23123253/221008849-91b5dd7c-2023-48a7-a347-c7ead54e863e.png)
Additional Information:
I often use
poisson_rng
in the transformed data block to implement a simulation of mock data from a Poisson process.Current Version:
I'm curious if others are able to reproduce this, or I am missing something. In case this is the wrong place to post the issue, I can also move it elsewhere. Thanks a lot!