# 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.