I have a question about Monte Carlo method. My Matlab model uses a parallel computing to accelerate the calculations with 32 cores.
The main goal is to obtain, after Monte Carlo processing, a Gaussian distribution, then to recover standard deviation. This distribution is provided from 100 000 draws in an uniform density.
Is that the same between:
to use 100 000 draws in an uniform density to obtain the output Gaussian distribution
and
to use 100 000 draws /32 cores in 32 uniform densities (the same that previous uniform density) to obtain 32 output single distributions and to add them to obtain the final Gaussian distribution? As showed in the picture.
Matlab gives the same result, but statistically is it right? Did you have any article that mentioned this similarly?
Thank you.
If you are using a pure Monte Carlo method then the answer is YES.
One of the conditions of the Monte Carlo processes is that all the events are independent, random and respond to the same probabilistic model. If this is your case, then your statement is totally right. Furthermore, You could run one single draw 100.000 times and adding them you should obtain the same result, as each draw is independent from the others.
Well, since the sum of two Gaussian distributions with the same mean and standard derivation is itself (up to a rescaling) a Gaussian distribution with the same mean and standard deviation, I don't see anything wrong with running the simulation on 32-cores (assuming I understood your problem description correctly).
What would concern me, however, is the independence of the random number generators. For instance, if the same RNG is used with the same seed value, instead of getting 32 independent sets of samples, you get the same set of samples repeated 32 times. Using a different seed but the same RNG may not be of much help either. First, the RNG may have a sequence length that is less than 100,000, so you may get the same set of numbers in a slightly different order. And even if the sequence length is not an issue, there is a possibility that the sequences generated on the different cores will be correlated.
One possibility would be to use a cryptographically strong RNG, but this may have an unacceptable cost on performance. Another thing I'd consider is the use of a custom, parameterized RNG with different sets of (well tested) parameters (not just seeds) for each core. Or perhaps even pre-generate a good set of 100,000 random numbers using a good RNG, and let each of the cores consume 3,125 of these.
I can only comment that the random generators are in fact pseudorandom generators and it a science on its own to use them correctly. Beginners often do the same mistake, for example, as mentioned above, running the same sequence of pseudorandom numbers or a shifted sequences on multiple cores or computers. Another mistake is done in generating random isotropic vectors by simply randomizing \phi uniformly in
Just a remark for testing your results: If you have a good analysis software with a nice graphical output and a non-linear fitting option (I use OriginLab), you can do a fit of your output data with a Gaussian function, knowing that the standard deviation of counts in every bin (n_i) should be \sqr(n_i). If your result is just a sum of identical, low-statistics, Gaussian data, the reduced \chi^2 will be far away from 1 (usually very high) and the analysis software would consider the result as inacceptable. This I do when getting data from instruments that do not seem to work properly but technicians are not sure. A high reduced \chi^2 of the Gaussian fit is a proof that something is wrong with the device and that the observed fluctuations do not correspond to the natural radioactive-decay fluctuations but may be caused by temporary interference by other devices nearby working with fixed-frequency signals (e.g. power converters). You can imagine how packed the instrumentation in the LHC is. :-)
You will need to initialize all the pseudo random number generators with independent high entropy seeds (a good practice regardless). If they all have the same seed then each will produce an identical set of values and so if that distribution is "not quite Gaussian" due to insufficient samples then the final distribution will simply be that "not quite Gaussian" distribution scaled by the number of cores.
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics
- Sampling Design
- Sample Size