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

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


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


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

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

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

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

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

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

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

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

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

I found the answer from the book Numerical Recipes in C, but unfortunately the code in that book is copyrighted. So, I re-derived a similar code that is slightly different but does the same thing. The idea is that we select a random number from distribution $g(x)$ that has the property that $f(x) < c g(x)$ always with a finite constant $c$, where $f(x)$ ...


1

So one thing you can do is implicitly represent a structured region, like a 3D rectangular volume, as a set of chopped up sub-volumes where each sub-volume automatically enforces the spatial constraints you require. You will do this implicitly by just storing how many pieces to slice up each dimension into individually and referencing each sub-volume using a ...


1

In the loop over your instances of 100 random values for $k$, you could write the value for $k$ into a file, and call your finite element code; the finite element code could then open the file, read the value for $k$ from it and then use it for whatever computation you may want to do with it. As for your second question, visualizing uncertainty is difficult....


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


1

Psudorandom numbers generated via methods like LCG are generated via a discrete dynamical system (i.e. a system where $x_n = f(x_{n-1})$). For these types of systems, tools such as Ergodic Theory can be applied to learn the the long-term distribution of such processes. This is one way to get the long-run probability distribution which generalizes well and ...


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