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

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 '17 at 18:24
• Sorry that I cannot remember its name. But I know this method is always correct. – user158565 May 24 '17 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. – Glen_b May 24 '17 at 22:38