13

To the best of my knowledge, Numpy does not support independent streams. Indeed, getting independent streams from the Mersenne Twister (Pythons RNG) is notoriously difficult although it can be done. Consider using the RandomGen package. It is fully compatible with Numpy, and provides you with the PCG64 generator, supporting up to $2^{63}$ independent ...


9

Like you say, using the Mersenne Twister for parallel computations is almost always done incorrectly, as the correct method is tricky to implement. By far the easiest and best answer would be to move away from the Mersenne Twister entirely, and use something like the PCG family, which provides multiple streams out of the box. The Mersenne Twister is known ...


7

A little playing with the sequence of numbers generated by the C code shows that the sequence is $z_{i+1}=5z_{i}+273 \mod 2^{16}$ This is a linear congruential generator (LCG). It's easy to show that this LCG has full period (See theorem 7.1 in Law's Simulation Modeling and Analysis, 5th ed. and check the three conditions.) I can't find the generator ...


6

To the best of my knowledge, no, but maybe other people here know the field more intimately. My knowledge comes primarily from developing Monte Carlo codes in physics. Knuth, in volume 2 of his Art of Computer Programming, states that Metropolis, using the middle-square method on 20 bits, found 13 cycles to which the method would always degenerate, the ...


6

I suspect the reason why generating the random numbers on the fly is slower for you is due to the rather large state of the Mersenne Twister. Switching to something like the PCG or XorShift+ random number generator would have several advantages for you: Higher quality of randomness (Mersenne Twister fails several tests for randomness) Smaller state, so ...


5

The initial velocities are drawn from a Gaussian distribution with variance $$\sigma_i^2=\frac{k_{\textrm{B}}T}{m_i},$$ where $k_{\textrm{B}}$ denotes Boltzmann's constant, $T$ is the temperature and $m_i$ is the mass of the $i^{\textrm{th}}$ particle. Thus, the problem boils down to generate random numbers from a gaussian distribution using uniformly ...


5

You definitely want to derandomize your program during development. Otherwise you will not be able to debug it since problems are not reproducible. At the same time, once you know the algorithm is working, you need to run it for multiple seeds or with different random number generators to ensure that your results (such as ensemble averages, standard ...


5

It is not a problem if one processor generates a number that has already appeared on a different processor. It would, however, be a problem if the two generated whole sequences that are similar. This is unlikely, if you start with different random seeds, since the number of numbers you get from a RNG before it starts to repeat itself (the cycle length) is so ...


5

A random number generator will give you random numbers that you can tweak to be between zero and 2^n and consequently it will allow you to sample random locations in your interval. These numbers may, in principle, repeat themselves -- so they are not a traversal, but that may not matter for the following reason: If $n$ is large (in the thousands), you will ...


4

I'd imagine that most users of random number generators are ultimately interested in floating-point values. This is why the Double precision SIMD-oriented Fast Mersenne Twister (dSFMT) exists. However, there is newer C code for the WELL RNG that returns unsigned long values. Looking at the code, it appears that the earlier version was casting unsigned long ...


4

As noted in the comments by Kirill, the y-axes of the two plots are very different. And if they are rescaled accordingly, the boxes will certainly look very similar, if not identical. Therefore, it is very reasonable to conclude that the raw simulation result in data.dat coming from your C++ code is correct, no matter what seed has been used for random-...


4

The random number generator is rarely the limiting factor in computational science. RNGs are usually quite simple and fast, a few dozen instructions, really. If you are doing anything even remotely complicated in your code with these random numbers, then the bottleneck is there.


3

If you want to use MT, you can use SFMT as your PRNG and SFMT jump to generate multiple streams. You can simply initialise MT with one seed, and then jump ahead by e.g. $1 \cdot 10^{60}$, $2 \cdot 10^{60}$, $3 \cdot 10^{60}$ … steps to generate multiple streams. Jumping is somewhat expensive, but you only need to do it once when you initialise your PRNGs.


3

Testing whether or not the mean is correct, or even if the histogram of your generated random variants "looks" like a certain distribution is not sufficient. Stick with much more rigorous test suites such as TestU01 or Diehard. Also, you really only have TWO random numbers in each row, because of the constraint that they sum to 1. This requires more ...


3

I'm not sure if these were added recently but it seems like there are now easy ways to generate random numbers quickly without too much overhead. From this article about Monte Carlo simulations in cython we can do from libc.stdlib cimport rand, RAND_MAX r = 1 + int(rand()/(RAND_MAX*6.0)) # random integer 1,...,6 As far as I understand you don't need to do ...


3

But these seem pretty old criticisms of the rand function This may be a nitpick, but I want to point out what I think is a flaw in this logic. Compilers are often extremely conservative about changing program behaviour, even when that behaviour (foolishly) depends on implementation details. This may or may not be true for the big compilers you’re familiar ...


2

If the data object you are trying to map is large enough, then just taking bit patterns is probably enough. For example, you might simply xor all the bytes of your object, giving you a number between 0 and 255 which you can then map onto the reals between 0 and 1.


2

The first thing that comes to my mind is the following: import random def randomy(obj): """ Map objects to uniformly distributed random numbers in [0, 1]. The input must be hash-able. You will get the same output for the same input during all calls in the same program run, but outputs may vary between different program runs. """ ...


2

This answer is late in coming, but you should have a look at SPRNG. It's specifically designed for scalability in parallel and supports a handful of types of PRNGs.


2

A simple idea for spreading a typical sequential RNG over a decent number of threads is to have a single thread advance the seed as fast as possible and send only every thousandth or so seed out to memory. Then have each of your other threads pick up one of these spaced out reference seeds and process the 1000 values in that block, i.e. regenerate again the ...


2

If you're interested in the distribution for fixed values of $a,b,p,q$, you can compute and tabulate the cumulative distribution by quadrature. Since it's only a 1d function for fixed parameters, you can compute the cumulative function for fixed values of $x$ with pretty high accuracy, connect these points by a polynomial approximation, and then use this ...


2

Most random numbers are in fact pseudo-random: they use a seed that is fed to a perfectly deterministic algorithm which can then generate series of numbers which are periodic but with a very large and complicated period so they usually appear random. This has drawbacks of course since it's not really random, but also advantages, among which is ...


2

For the sake of your example, let's say RAND_MAX=12 and i=17. Then do the following procedure: Choose two random numbers $r_1,r_2$ and combine them to a single random number uniformly distributed in $[0,144)$ by computing $r=r_1*12+r_2$. This is of course the wrong interval. You get a uniformly distributed random number in $[0,17)$ by repeating this process ...


2

Instead of fixing the provided generator use a reasonable modern choice. It will be faster and have better statistical quality. Possible examples include the various variants of xorshift+, xorshift* and PCG. To directly respond to the question asked. You can generate a sample, mask out the maximal power of two bits available, generate another and shift ...


2

To permute a list, there is no need to enumerate the n! possibilities. Shuffling the list should be sufficient. Pseudo code (from wiki article) for shuffling would be for i from n−1 downto 1 do j ← random integer such that 0 ≤ j ≤ i exchange a[j] and a[i] For c++11, std::shuffle along with a good RNG should work. This wiki article was my source. ...


2

You can do some spot checking. Testing for all $n!$ possibilities is impractical, but you can make some unit tests to make sure that after doing sufficiently many shuffles certain properties of the distribution are correct. For example, make the array $[1,2,\ldots,n]$ and check (for $x=2,5,7,12$) Is the mean value for the index of $x$ after shuffling ...


2

If you can use C++11 you can use the built-in Poisson distribution #include <iostream> #include <random> static std::random_device rd; static std::mt19937 gen(rd()); int main() { std::poisson_distribution<int> pd(5); for (int i = 0; i < 10000; ++i) std::cout << pd(gen) << '\n'; } I mean, come on, it doesn't get ...


2

Welcome to SE SciComp. First of all, I would suggest using Jupyter so that you have access to IPython and its nice timing magics (see %time and %timeit magic function). These magics take into account that you need to run the code a couple of times to get reliable measurements. If you need even more in-depth comparison, you also need to take into account ...


1

This exact problem (not generating the same random number sequences) is described in John D. Cook's blog post Random number generator seed mistakes. The solution is simply to ensure that processes are started with an unique seed. Which you also mention as a possibility in your question. I think your figure might be a bit misleading as most random number ...


1

Really only you can answer the question about simulation bias and if it is acceptable in your application. The standard procedure I use: Set a pseudo random sequence as a benchmark (standard Monte Carlo) using a high # of simulations (in risk management 10,000 is often used, in other fields 100,000 to 1M may be used). Run your RNG over the same input data ...


Only top voted, non community-wiki answers of a minimum length are eligible