2
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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 three trapezoid methods
    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.]

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks, this seems a form of rejection sampling, right? $\endgroup$ – rano May 24 '17 at 18:24
  • $\begingroup$ Sorry that I cannot remember its name. But I know this method is always correct. $\endgroup$ – user158565 May 24 '17 at 18:46
  • 1
    $\begingroup$ @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. $\endgroup$ – Glen_b May 24 '17 at 22:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.