Sampling from an Arbitrary Density Math, CS, Data
One way to sample from a known probability density function (pdf) is to use inverse transform sampling. First, you integrate the pdf to get the cumulative distribution function (cdf). Next, you find the inverse of the cdf. Finally, apply this inverse cdf to each number in a sample of Uniform(0,1) observations.
Suppose we wish to sample from the pdf
The cdf is
For $u \in [0,1]$, its inverse is
Finally, we can apply the technique. Below, I generate 100 standard uniform observations and transform each using the inverse cdf.
R Function for an Arbitrary pdf of One Variable
The R code below uses some of R’s built-in numerical methods to accomplish the inverse transform sampling technique for any arbitrary pdf that it is given. I’m sure there are some ugly pdfs for which this function wouldn’t work, but it works fine for typical densities.
We can use
samplepdf to sample from $h$, which we defined earlier.
Computational Speed of the Algorithm
samplepdf is convenient, it is also computationally slow. The algorithm is $O(n)$, so it is easy to estimate the time required for a large sample size. Simply fit a linear model to the times required for a set of small sample sizes, as shown below.
At the beginning of this document, we found the inverse cdf analytically. The R function that used this analytical result can generate a sample of size 1000 in the blink of an eye, while the
samplepdf function took about two and a half minutes on my slow netbook. However, if you give the lower and upper endpoints as in
samplepdf(1000, h, m = .5, spdf.lower = .5, spdf.upper = 1), the function only takes about ten seconds to run.
Another option for sampling that doesn’t require finding the inverse cdf is rejection sampling.
blog comments powered by Disqus