How to sample numerically from an arbitrary smooth distribution?

Question

I'm given a smooth probability density function via its values on a reasonable fine grid. I assume that cubic spline interpolation (or cubic spline interpolation of the logarithm of the density) will be sufficient to evaluate it at arbitrary points with high accuracy. I wonder how to generate random numbers that reproduce this distribution.

My first shot was to approximate the cumulative distribution function of this distribution by a piecewise linear function $F$ (on the original grid), draw a number $r$ from $[0,1)$ uniformly at random, and take the $x$ with $F(x)=r$. However, I noticed that the accuracy of my final results is not great, and I suspect that I lose accuracy because the piecewise constant probability density of my numerical random variable doesn't approximate the real smooth probability density function well enough. What options do I have?

Here are some of my ideas:

1. Go to the library and look for a book about Monte Carlo simulation. Or try to ask an expert.
2. Integrate the cubic spline analytically, which gives a piecewise quartic function $F$. There would still be an analytic formula for the $x$ with $F(x)=r$, but it will probably be complicated to implement and slow to evaluate.
3. Approximate the smooth probability density function by a piecewise linear function, which gives a piecewise quadratic function $F$. The analytic formula for the $x$ with $F(x)=r$ should be simple to implement and reasonably fast to evaluate.
4. Approximate the logarithm of the smooth probability density function by a piecewise linear function, which gives a piecewise "simple" analytic function $F$. The analytic formula for the $x$ with $F(x)=r$ should be simple to implement and reasonably fast to evaluate.
5. Approximate the smooth probability density function $g$ by a piecewise constant function $f$ such that $g \leq 1.1 f$. Now use rejection sampling by first sampling $x$ via $F(x)=r_1$, and then rejecting $x$ if $g(x) < 1.1 f(x) r_2$.
6. Approximate $F^{-1}(r)$ by a suitable piecewise analytic function. But what does suitable mean here?

Answer 1

If your PDF is bounded, you could try approximating its inverse with a high-degree polynomial interpolant. This is usually considered a bad thing, but that's just a myth.

Some things to keep in mind:

• Instead of using an equispaced grid, interpolate at the Chebyshev nodes of the first kind, i.e. $x_i = \cos\left(\pi\frac{2i-1}{2N}\right)$, for $F(x)$ defined in $[-\infty,\infty]$, or second kind, i.e. $x_i = \cos\left(\pi\frac{i-1}{N-1}\right)$, for $F(x)$ defined on a finite interval.

• If $F(x)$ is on an infinite interval, don't interpolate $F^{-1}(r)$, as it will have singularities at the endpoints, but interpolate $F^{-1}(r)/(r^2-r)$, as this will cancel-out the singularities at $r=0$ and $r=1$. Using the Chebyshev nodes of the first kind will avoid evaluating $F^{-1}(r)$ at these singular points.

• You can evaluate your interpolant using Barycentric interpolation. Note that if you evaluated $F^{-1}(r)$ on Chebyshev nodes of the first or second kind, the Barycentric weights $w_j$ have closed-form expressions.

• For a faster evaluation that vectorizes well, use a Vandermonde-like matrix $V$ with $V_{ij}=T_{j}(x_i)$ to compute the Chebyshev coefficients of your interpolant once (if you used the Chebyshev nodes, $V$ should be well conditioned) and use Clenshaw's algorithm to evaluate it for more than one $r$ at a time.

The method described here is more or less what the Chebfun system does (disclaimer: I used to be part of the Chebfun developer team). Most of the basic Chebyshev technology is described in Nick Trefethen's book "Approximation Theory and Approximation Practice", of which the first six chapters are available online.

Answer 2

I would go with your option 3.

That said, it would have helped if you elaborated on your statement "However, I noticed that the accuracy of my final results is not great". I say so because if the mesh on which your PDF is defined is fine enough, then I see no reason why your approach should not work. What I would do is try to debug things by starting with a PDF you know analytically, say a Gaussian, and evaluate the steps you do one by one. For example, start with a very fine mesh and piecewise constant approximation -- does the resulting set of samples look ok? If not, does it get better by using a piecewise linear approximation? If not, then the error must be somewhere else. Etc.

Answer 3

I have solved my problem now. The reason why I lost accuracy doesn't even occur in the question. Let's first address the proposed solution ideas from the question:

1. Trying to learn more was a good idea in this case.
2. Analytical formulas are attractive when they are reasonably simple. This is not the case here.
3. Straightforward linear interpolation sounds like a good idea, at least for verification. Perhaps I will implement it one day, it can't be too difficult. See also the answer by Wolfgang Bangerth.
4. Note sure. I have done something related instead. I approximate the smooth probability density function $f(x)$ by a piecewise analytical function of the form $a \cdot x^b$. Its antiderivative $\frac{a}{b+1}x^{b+1}$ has the same simple form, and the moments $\int f(x) x^n dx$ also lead to integrals of the same form.
5. Rejection methods shouldn't be outright rejected, but I haven't tried this one.
6. Approximating $F^{-1}(r)$ is one of the "correct" proceedings. What does suitable mean here? The inverse cumulative distribution function is monotonically increasing, and a near zero probability density over an extended range translates into a very steep slope of $F^{-1}(r)$. Nonuniform rational interpolation is one option able to cope with these features. The nonuniform grid leads to a $O(\log n)$ effort for a single function evaluation, but this is normally still fast enough in practice. See also the answer by Pedro. (However, the $O(n)$ effort for a single function evaluation with Chebyshev interpolation instead of a $O(1)$ or $O(\log n)$ effort has made me uneasy in the past.)

---

But what about the implicit "real" question?

However, I noticed that the accuracy of my final results is not great, and I suspect that I lose accuracy because ...

I actually wasn't just given a single probability density function, but a one parameter family of probability density functions, tabulated on a sufficiently fine grid relative to the one parameter. I handled this by linear interpolation between the inverse cumulative distribution functions I had precomputed. However, even in case the probability density function is well approximated by a linear interpolation between the tabulated probability density functions, my above proceeding can lead to unacceptable errors.

The "correct" solution is much simpler. It uses the composition method. First draw a number uniform at random to select between the available precomputed distributions, and then sample the value from the selected distribution. This actually generates the "exactly correct" distribution, in case the probability density function is really given as a weighted sum of available precomputed distributions.