7

I am investigating whether it is possible to have two sets of software agree on a sequence of produced pseudo-random numbers. I am as interested in understanding all the possible points of divergence as I am in actually finding a way to get them to agree.

Why? I work in a data shop that uses many different software packages (Stata, R, Python, SAS, probably others). There has recently been interest in QCing outputs by replicating processes in another language. For any process that involves random numbers, it would be helpful if we could provide a series of steps ("set this option", etc.) that allow the two packages to agree. If that's not feasible, I'd like to be able to articulate where are the failure points.

A simple example:

Both R and Python's default random number generator is Mersenne-Twister. I set them to the same seed and try to sample from and also look at the "state" of the PRNG. Neither value agrees.

R (3.2.3, 64-bit):

set.seed(20160201)
.Random.seed
sample(c(1, 2, 3, 4, 5))

Python (3.5.1, 64-bit):

import random

random.seed(20160201)
random.getstate()
random.sample([1, 2, 3, 4, 5], 5)
ebyerly
  • 662
  • 5
  • 14
  • 4
    Why not implement a standard algorithm yourself? I know at least for Python, there's no guarantee the RNG implementation won't change. It *has* changed, actually; it used to use Wichmann-Hill, and a number of methods used to rely on `random.random()` when they really shouldn't have. (Python 2.7 still does misuse `random.random()`, actually, which can cause some [weird issues](http://stackoverflow.com/questions/33602014/random-number-in-the-range-1-to-sys-maxsize-is-always-1-mod-210).) – user2357112 Feb 01 '16 at 21:03
  • 2
    I was about to suggest the same -- you could just port the simple + fast [RcppZiggurat](https://cran.rstudio.com/web/packages/RcppZiggurat/index.html) to Python. – Dirk Eddelbuettel Feb 01 '16 at 21:08

1 Answers1

1

Old question, but maybe useful to some future reader: As alluded in the comments, your best bet is to implement this your self and provide interfaces for the different environments such that for a given seed the same results are returned. Why is that necessary? You used "sampling" as an example. There are several steps involved.

  1. Seeding is a non-trivial process. For example R goes as far as to further scramble the provided seed. So unless you user tools use the same method, they will end up with a different seed even when the user supplies the same value.

  2. The actual RNG: Even though in both cases Mersenne-Twister might be used, is it really the same version that is used? R uses a 32bit MT. Maybe Python uses a 64bit version?

  3. Most RNGs give you an unsigned integer (nowadays typically 32 or 64bits). But you will need some distribution of random numbers, e.g. for sampling you need random integers within a given range. There are many methods to go from the integers produced by the RNG to those needed for sampling. In the case of R, you do not even have access to the output value of the RNG. The most fundamental function is R_unif which returns a double in [0, 1). Again, how to generate such a double is not universally agreed on. And if you need other distribution functions (normal, exponential, ...) you will find quite a few different algorithms for them.

Overall there are to many places where (subtle) differences can creep in.

Ralf Stubner
  • 26,263
  • 3
  • 40
  • 75