low-discrepancy sequences for sampling of distributions

Question

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

Answer 1

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.

Answer 2

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.