# Generate a Monte Carlo sample from a PDF defined by a Fourier Series

I have a probability distribution (PDF) defined by a Fourier series.. actually it's a purely cosine series over a known range. The PDF quite smooth, so most of the power is in the low 5 or so frequencies. The function is (like a PDF must be) always positive. It's well defined everywhere.

I only need a single sample, so tabulation or function fitting is likely too expensive; that'd pay off if I needed thousands of samples, but I need just one, then I'll throw the PDF away and generate a new one. (This is all part of a Markov chain walk, and the transitions are defined by Fourier coefficients).

So, I want to efficiently generate one sample from this PDF. What's the best strategy?

One method that works but is slow is to use a von Neuman rejection method. The PDF is bounded from above by the sum of the absolute values of all the sine coefficients. I can generate a random point x in U(0,1), evaluate the series for that x, then accept the point with a probability based on the ratio of f(x) and the upper bound. This will work but it's super-slow!

The second strategy would be to generate a U(0,1) sample and feed it through the inverse CDF, but that involves solving a nonlinear equation for x, which is iterative and even slower.

Are there any other strategies I can experiment with?

-
You speak of "evaluating the series" -- it seems that what you have is actually a Fourier sum? – joriki Jun 23 '11 at 7:14

You could divide the interval into subintervals, for each subinterval precalculate the minimal and maximal values of each cosine there, and then when you have a distribution you can use that to calculate better bounds in each interval. Then you need to sum the areas under all the bounds, generate one random number to choose a subinterval according to the summed areas, and then apply the rejection method within the subinterval. However, if most of the power is in the lowest $5$ modes, I doubt that this will lead to a significant improvement, if any; in fact I doubt that in that case you'll be able to significantly improve on the simple rejection method with any method.

-

The solution has a nice elegance about it.

Since this is a Fourier series, there's both cosine and sine terms over the range -pi/2 to pi/2. Break the sum up into the two halves. Notice the cosine terms are symmetric around 0, the sine terms are antisymmetric. So if we can first properly sample a point x from the cosine series alone from 0 to pi/2, we can use the simple evaluation of F(x) versus F(-x), using the ratio to decide which point to accept. The sine terms basically let you choose between the two alternatives.

But this still leaves the problem of sampling from a (simpler) cosine series from 0 to pi/2. Well, we can transform that recursively into symmetric and antisymmetric terms over the range 0 to pi/2! You can do this by shifting and stretching the range, making it a Fourier series again (which can now be sampled recursively), or by just keeping track of odd and even cosine terms... the methods end up being equivalent.

So each recursive level breaks down the sampling to a self-similar problem but with half the coefficients. Since Fourier coefficients rapidly decrease (assuming the PDF is pretty smooth) you'll likely end up with the simple problem of sampling from a single cosine lobe with ignorably small noise added to the lobe. Sampling from a cosine is easy (asin(x)), thus terminating the recursion.

In practice, this all works stably and fast, especially if you have just 4 or 8 terms.

-