Sampling from a piecewise linear probability density function

Question

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?

Answer 1

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

Answer 2

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.