2

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?

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
Fortranner
  • 2,525
  • 2
  • 22
  • 25
  • 2
    Google suggested [this](http://stackoverflow.com/q/1594121/324364). – joran May 09 '14 at 17:50
  • Thanks, that question and mine are indeed the same. – Fortranner May 09 '14 at 19:11
  • 1
    See [this](http://stackoverflow.com/questions/20508400/generating-random-sample-from-the-quantiles-of-unknown-density-in-r/20509483#20509483) and [this](http://stackoverflow.com/questions/23258482/use-inverse-cdf-to-generate-random-variable-in-r/23259123#23259123), and [this](http://stackoverflow.com/questions/23223548/how-can-i-generate-a-sample-from-a-log-normal-distribution-with-pareto-tail-in-r/23228354#23228354). – jlhoward May 09 '14 at 19:27
  • Actually, your question and the Google reference in @joran's comment *not* the same. Your are asking if there is a package in R to do this. There was no package in 2009, when that question was asked and answered. But there is now (see my response below). IMO this is a *perfect example* of why embargoing questions as "duplicated" is a really horrible practice (although no one is suggesting that for your question). – jlhoward May 09 '14 at 19:49
  • @jlhoward , good point, although then we get into a "close because asking to recommend a package" situation ... (it's not entirely clear to me that this close criterion should really be as applicable for [r] as it is for other tags ...) – Ben Bolker May 09 '14 at 19:50
  • @jlhoward In my defense, I simply provided what assistance I could in a comment. I never voted to close until the OP _explicitly_ said that the linked question was a duplicate. Frankly, the best solution would be for your answer to be on the old question, and for this to be closed as a duplicate of that one. – joran May 09 '14 at 19:52
  • @BenBolker Yes I know. Still I think expecting someone to come across the `distr` package by accident is a pretty big lift. If you google "distr in R" you get nothing about this package. – jlhoward May 09 '14 at 19:52
  • @joran Sorry if this came across as directed at you personally - it was definitely *not* intended that way. I just believe embargoing questions as duplicates is a really bad practice in general, and this was (IMO) one good example of why. – jlhoward May 09 '14 at 19:56

2 Answers2

7

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.

jlhoward
  • 58,004
  • 7
  • 97
  • 140
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")

enter image description here

josliber
  • 43,891
  • 12
  • 98
  • 133