# low-discrepancy sequences for sampling of distributions

Hi everybody, I am writing of physics simulation that traces charged particles. I need to set up these particles in phase space ( 3 space dimension + 3 momentum dimensions) following distributions. For example I need a Gaussian distribution in the xy plane and a Gaussian in the z direction and same in momentum space. So i need to sample the 6 dimensional parameter space effectively to reduce the number of particles. I found low-discrepancy sequences and think they allow a $\frac{1}{N}$ scaling instead of the $\frac{1}{\sqrt(N)}$ for random sampling.

I am not sure how to apply them. I thought of creating a 6 dimensional Halton sequence and then plugging them into the cumulative distribution function.

I tried a simpler example to create only a Gaussian distribution in the xy plane using a 2d Halton sequence as input for the Box–Muller transform to create a Gaussian distribution. But the histograms for the x and y direction to not look as evenly distributed as I would like them to be. Am I using the concept of low-discrepancy sequences right?

Thanks for any help

• Did you try physics.stackexchange.com, stats.stackexchange.com, or stackoverflow.com? – Harry Gindi Feb 9 '11 at 11:02
• @Harry, I don't know an enormous amount about this topic, but I think I know enough to say that it's genuine mathematics, certainly not physics, on the mathy side of statistics. – Gerry Myerson Feb 9 '11 at 12:00
• i didn't try the other sites, but I think it is a rather mathematical question so I think this is the right place. – Radfahrer Feb 9 '11 at 12:33
• I would recommend using Sobol' sequences as they tend to give better results in my experience. Second, using low-discrepancy sequences you do not want to use methods like Box-Muller to generate normals from the uniforms since, as pointed out by Tim, the points in the low-discrepancy sequences are not independent. Simply invert the cumulative normal density function (the Ziggurat method is an efficient method for this). – Apollo Feb 23 '11 at 15:56

Using Box Muller to generate normals from a uniform quasi random sequence indeed will give rather ragged marginal distributions that look similar to the distributions that a pseudorandom uniform distribution with Box Muller will generate. Much better is to use an explicit inverse normal cdf applied to the quasi random sequence. This guarantees that the empirical cumulative distribution function of the samples is a low discrepancy sequenc, and the marginals will now look like almost 'perfect' bell shaped normal distributions. There is no such guarantee after you apply the Box Muller transformation.

Let's see if I understand your question:

From your introduction I thought you need to sample 6 coordinates for a fixed number of particles from a prescribed (6-dim, i.e. multivariate) probability distribution as an initial condition (for a simulation of the evolution of a n-particle system following classical mechanics and electrodynamics). Therefore I don't understand this statement:

• "So i need to sample the 6 dimensional parameter space effectively to reduce the number of particles."

How does the number of particles depend on your sampling algorithm? Or are we talking about a stochastic evolution? (Maybe we are conflating "particle" in the sense of a point mass in physics and "particle" in the sense of a stochastic particle filter?)

• "I found low-discrepancy sequences and think they allow a $\frac{1}{N}$ scaling instead of the $\frac{1}{\sqrt{N}}$ for random sampling."

The latter is the order of the error of a Monte-Carlo integration, so now I suspect I misinterpreted your first statements, did I?

• "I thought of creating a 6 dimensional Halton sequence and then plugging them into the cumulative distribution function."

6 dimensional Halton sequence:

The usual strategy to sample from a prescribed distribution is to sample from a uniform distribution on $[0, 1]$ and to apply an appropriate transformation. If you use a deterministic algorithm to generate samples from a distribution, you'll always have the problem that the generated numbers are correlated. It does not matter, from a theoretical viewpoint, if you initialize your algorithm once with a certain seed and use it for the six coordinates of your first particle in sequence, then for the second etc. or if you initialize six versions of your algorithm with six different seeds, and use one version to generate one coordinate of your first particle etc.

Is there a specific reason why you chose Halton sequences instead of some random number generators from any of the standard introductory sources like "numerical recipes"?

From a theoretical viewpoint, both algorithms don't generate truly random numbers. From a practical viewpoint, both algorithm may or may not produce artefacts in your simulation, you'll have to evaluate this from case to case, but there is no reason to believe a priori that you fare better with the concurrent use of six versions of the algorithm.

• "But the histograms for the x and y direction to not look as evenly distributed as I would like them to be."

Since no deterministic algorithm produces truly random numbers, the best one can say about an algorithm is if it passes a fixed set of statistical tests. One of the first tests everybody would use is a test that tests the plausibility of histograms, a chi square test. If you implement a well known pseudo-random generator that fails such a test, it is rather safe to say that your implementation has a bug.

For an up to date source on random number generation, including multivariate normal distributions, see

• James E. Gentle: "Random Number Generation and Monte Carlo Methods"

For more details and references, see random number generator on the Azimuth project.

• sorry I didn't formulate the question right. I don't have a single 6 dimensional distribution function. What I want is to have a product of different distribution functions. So for example a Gauss distribution in xy plane and a linear distribution in the z direction. So i need to create different distributions. The particles are advanced by classical mechanics. I could use plain random numbers, but I would need to simulate more particles, to reduce the error. So if I understand it correctly I can use low-discrepancy sequences to sample my different distributions. – Radfahrer Feb 9 '11 at 15:15
• i am trying to archive something similar like this: google.de/… in chapter 2.6 Initial particle distribution (page 48). I know the use the Hammersley sequence, but don't know how – Radfahrer Feb 9 '11 at 15:26
• @Radfahrer: A particle has 6 coordinates which are random variables and have a joint distribution function, ergo you have a 6-dim-distribution function - what you are saying is that some of these random variables are independent and the 6-dim distribution factors into the product of a two dim in xy coordinates and a one-dim in the z-coordinate etc. My main problem now is to understand, when you say "I would need to simulate more particles, to reduce the error", the error of what? I'll look at the reference you linked to... – Tim van Beek Feb 10 '11 at 8:12
• @Radfahrer: THe reference you linked to is about a simulation of macroparticles in prescribed external electromagnetic fields and boundary conditions, in order to run simulations of new accelerator designs - you'll always have to simulate a certain amount of macroparticles to get representative results, mainly depending on the projected beams intensity - no matter what kind of random number generator you use to generate your initial distribution. Just try one from the standard sources, like "numerical recipes". – Tim van Beek Feb 10 '11 at 8:39
• I want to calculate the beam width as standard deviation of the particles. This of course depends on the initial distributions. To get good results I need more particles. In the link they have those nice histograms for only 200 particles that represent the Gaussian distribution really good, and a plot with random numbers, which is not very good. This nice histogram is created with a low-discrepancy sequence. I simulate single electron beams. I just use more particles for sampling the initial distributions to get meaningful statistics for the beam width. – Radfahrer Feb 10 '11 at 9:09