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 ...


7

Cython makes code faster by removing the type ambiguity. Since random.py is a pure python module, you can just copy it and add the types to the functions you need. Then cython can optimize the dynamic overhead away.


6

If I understood your question correctly: Subdivide the interval $[0,1]$ in $N$ segments, each having a width proportional to your probabilities $p_i$ (where $\sum_{i=1}^N p_i = 1$). Use a pseudo-random-number-generator using the uniform probability density distribution on $[0,1]$ to generate the number $x$. Determine in which segment $x$ is located; this ...


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

Following aterrel's suggestion, you could use pyximport to automatically compile the random module: import pyximport pyximport.install(pyimport=True) import random However, this still will not make it as fast as it would be if you declared static types for the variables in Cython.


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-...


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 think the book Numerical Recipes--The Art of Scientific Computing,William H et al talks about each PRNG's strengths and weaknesses in chapter 7. If you want to test different PRNGs by yourself, you can try the TestU01 suite, which is developed by L'Ecuyer and Richard Simard, this suite contains 3 predefined batteries(SmallCrush, Crush, and BigCrush), ...


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

I wanted to chime in to say that, I've recently gone through this process with a simulation and I should note that using Fortuna is not out of the question if it is really necessary. In our case, we were concerned that MT's entropy wasn't high enough which would translate in our simulation to a bias. So for our simulation we used Fortuna pulling about 65 ...


2

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 ...


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

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

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

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

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 ...


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