How to sample numerically from an arbitrary smooth distribution?

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?
• Could you give a few more details as to your environment? Are you doing this in compiled code, or in some system such as Matlab or Python? Sep 30 '13 at 9:37
• @Pedro This is compiled code, more precisely C++. A complete simulation takes some minutes on a modern Intel CPU with 12 threads, during which several million values are drawn from such distributions. The preprocessing time is currently completely neglegible compared to the time taken by std::lower_bound for drawing $x$ via $F(x)=r$. Sep 30 '13 at 10:50

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.

• Pretty sure you wanted to type F^{−1}(r)*(r_2−r) instead of F^{−1}(r)/(r_2−r). Otherwise you add "more" singularities. Or am I having a brain fart?
– oli
Aug 21 '20 at 15:55
• And although Chebyshev interpolants are in general very useful, in this case I personally prefer a low order equidistant spline. Usually one samples a lot of numbers, and such splines are faster in the evaluation with the usually negligible trade off of increased number of interpolation points, i.e. slower in the construction and more memory. If F^{-1} can be interpolated directly (when it has no singularities) one can chose a monotonic cubic spline to ensure monotonicity of the interpolant as well, which sometimes is useful.
– oli
Aug 22 '20 at 19:56

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.

• A added that I used the original grid for the piecewise linear approximation of the cumulative distribution function. However, this grid already has more than 600 points, and the final results also look more or less OK, just not perfect. I do have (various) approximate analytical expressions which I can use instead of the tabulated density, so I have quite some options for debugging things and evaluating the accuracy of my final results. My current feeling is that learning about commonly used sampling sampling schemes and using a better converging one is a good next step. Sep 30 '13 at 8:36
• 600 points may or may not be sufficient. It all depends on how peaked your distribution is -- for example, 600 samples between -100 and +100 is probably not enough if you are sampling the Gaussian distribution. Oct 1 '13 at 2:12
• The logarithm of the density is quite smooth, the density itself is "a bit" peaked at zero (the distribution "starts" at zero). Of course you may be right that other inconsistencies might be responsible for the loss of accuracy, but debugging randomized algorithms like Monte Carlo is always slightly challenging for me. I just wanted to learn about the area where I was lacking most knowledge, which I have done now. Oct 1 '13 at 22:39
• Fair enough. And yes, debugging MC is difficult. But if you know the distribution and it's simple enough, then it's not impossible within an hour or two to get a billion samples and reconstruct most quantities you may be interested in to pretty good accuracy. This helps find most systematic errors that may not be visible using just a few thousand samples due to the randomness of samples. Oct 2 '13 at 12:52

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.

• Keep in mind that an interpolation polynomial over $n$ points will have degree $n-1$, and thus you will probably need a much smaller $n$ than in the piecewise linear case. Nov 1 '13 at 10:58
• @Pedro What made me uneasy in the past was the case where the size of the domain determined the required number of Chebyshev interpolation points. The number of Chebyshev points was determined automatically, but the size of the domain depended on the problem and could become quite huge. Perhaps this observation suggests one additional possible remedy (in addition to the ones suggested in the linked question/answer): Keep the number of Chebyshev interpolation points bounded by some constant ($\approx 16$), and subdivide the domain in case more interpolation points are needed. What do you think? Nov 1 '13 at 11:48
• The number of points is not determined by the size of the domain, but by its complexity, i.e. by how difficult it is to approximate it by a polynomial. Recursive bisection is, in any case, a good option. Nov 1 '13 at 12:08
• @Pedro The domain is always an interval, there is no complexity here. My past case was optics related, and you roughly were required to have at least two points per wavelength. I used an analytical formula which determined exactly how many points I needed, but the two points per wavelength are enough to understand the problem. As soon as the size of the domain was many wavelengths (say $\approx 20$ wavelength), the speed of Chebychev interpolation became a potential issue. Nov 1 '13 at 12:17