Simulate from an (arbitrary) continuous probabilit

Question

This question already has an answer here:

• How do I best simulate an arbitrary univariate random variate using its probability function? 3 answers

For a normalized probability density function defined on the real line, for example

p(x) = (2/pi) * (1/(exp(x)+exp(-x))

(this is just an example; the solution should apply for any continuous PDF we can define) is there a package in R to simulate from the distribution? I am aware of R's built-in simulators for many distributions.

I could numerically compute the inverse cumulative distribution function at a set of quantiles, store them in a table, and use the table to map from uniform variates to variates from the desired distribution. Is there already a package that does this?

Answer 1

Here is a way using the distr package, which is designed for this.

library(distr)
p    <- function(x) (2/pi) * (1/(exp(x)+exp(-x)))  # probability density function
dist <-AbscontDistribution(d=p)  # signature for a dist with pdf ~ p
rdist <- r(dist)                 # function to create random variates from p

set.seed(1)                      # for reproduceable example
X <- rdist(1000)                 # sample from X ~ p
x <- seq(-10,10, .01)
hist(X, freq=F, breaks=50, xlim=c(-5,5))
lines(x,p(x),lty=2, col="red")

You can of course also do this is base R using the methodology described in any one of the links in the comments.

Answer 2

If this is the function that you're dealing with, you could just take the integral (or, if you're rusty on your integration rules like me, you could use a tool like Wolfram Alpha to do it for you).

In the case of the function provided, you can simulate with:

draw.val <- function(numdraw) log(tan(pi*runif(numdraw)/2))

A histogram confirms that we're sampling correctly:

hist(draw.val(10000), breaks=100, probability=T)
x <- seq(-10, 10, .001)
lines(x, (2/pi) * (1/(exp(x)+exp(-x))), col="red")