## 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 2011 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 2011 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 2011 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 2011 at 15:56

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.

-