# Sampling from a piecewise linear probability density function

I have a pdf modelled as a piecewise linear function that I can represent as the sequence of pairs of real numbers standing for the points in the of the piecewise line: $\{(x_{1}, y_{1}), \dots, (x_{n}, y_{n})\}$ such that each $x_{i}, x_{i+1}, i=1\dots n-1$ is a bin; the first two points are on the x axis, i.e. $y_{1}=y_{n}=0$, and the (trapezoid) area defined by the line is equal to 1 (it is a proper probability distribution function).

How can I procedurally sample from such a distribution? I feel like I can first determine one bin, e.g. $x_{s}, x_{s+1}$, by randomly selecting one with probability proportional to the trapezoid mass inside the bin, then how can I sample from the line defined by the two points $(x_{s}, y_{s}), ( x_{s+1}, y_{s+1})$?

Or shall I go through the cumulative distribution function? How can I derive it from the piecewise density representation of a sequence of points?

I believe your suggested method is likely to be the best way to do it in general.

step 1: sample an interval to generate from, using the discrete probability distribution of the trapezoid areas.

step 2: sample from the conditional distribution given you're within the trapezoid (i.e. as if it was scaled to be its own density).

There are a number of efficient algorithms for step 1. See: How to sample from a discrete distribution? or How to generate numbers based on an arbitrary discrete distribution?

Now step 2 can be done in a variety of ways.

It's probably easiest to just write a generic trapezoid sampler (though if you're writing within particular platforms you may already have one). e.g. assume you have a trapezoid over $$(0,1)$$ (with only one parameter, $$h$$ - the height at $$0$$, where $$0\leq h\leq 2$$), then scale and shift the result as needed for each segment. Trapezoids can be generated in a variety of ways, for example you might:

• use inverse cdf -- the cdf is quadratic, so not hard to invert; $$X = \frac{\sqrt{h^2\, +\, 4\,U\,(1-h)}\,-\,h}{2\,(1-h)}$$, but you need to check for $$h=1$$ and return $$U$$ in that case (and if h is likely to be very very close to 1 you may be better to rewrite the function so that it doesn't suffer from catastrophic cancellation).

• treat it as a mixture of a uniform plus a triangular; the triangular is straightforward in any of several ways.

• you can treat it as a mixture of two triangular distributions

• you could even use simple rejection sampling on each segment (in many practical cases of this problem it will be fairly efficient and it can never be worse than 50% rejection) Illustration of the last three trapezoid methods

If you needed many draws from this distribution, you'd precompute the information needed for the discrete sampler and precalculate all the $$h$$ and x-scaling values for each segment -- then generation from the distribution as a whole should proceed quickly.

There are numerous other ways to generate from a trapezoid not mentioned here. [Note that this piecewise trapezoid needn't be everywhere continuous; if you keep track of both ends for every segment the same approach would work just fine.]

• If you wanted to sample from some distribution could you use a piecewise linear approximation to the true density as a proposal distribution and then use the techniques in your post to sample from the proposal distribution (then the proposal distribution samples are used by some sampling technique to sample from the true distribution)? I was reading about sampling techniques from some lecture notes (based on PRML) and they never explained how you get the proposal distribution for the sampling methods, actually when I read PRML I don't recall them explaining where it comes from either. May 25 at 4:21

One inefficient but correct method:

Suppose $x_1 < x_2 <...<x_n$.

1: Generate two random numbers: $X$ ~ uniform($x_1,x_n$), $Y$ ~ uniform($0,\max(y_1,..., y_n)$).

2: Check if $Y<f(X)$, where $f$ is piece-wise linear function used to define your density. If yes, keep $X$, otherwise, discard X.

3: Repeat 1 and 2 until you achieve your sample size.

• Thanks, this seems a form of rejection sampling, right?
– rano
May 24, 2017 at 18:24
• Sorry that I cannot remember its name. But I know this method is always correct. May 24, 2017 at 18:46
• @rano yes, that's basic rejection sampling. In some cases it may be okay, but if $\max(y_1,\ldots,y_n)$ is a lot larger than say $E(f(X))=\int_{-\infty}^\infty f^2(x) \, dx$ its efficiency might be extremely low. May 24, 2017 at 22:38