Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Generate random variables from a distribution function using inverse sampling

I have a specific density function and I want to generate random variables knowing the expression of the density function.

For example, the density function is :

df=function(x) { - ((-a1/a2)*exp((x-a3)/a2))/(1+exp((x-a3)/a2))^2 }

From this expression I want to generate 1000 random elements with the same distribution.

I know I should use the inverse sampling method. For this, I use the CDF function of my PDF which is calculated as follows:

cdf=function(x) { 1 - a1/(1+exp((x-a3)/a2))

The idea is to generate uniformly distributed samples and then map them with my CDF functions to get an inverse mapping. Something like this:

random.generator<-function(n) sapply(runif(n),cdf) 

and then call it with the desired number of random variables to generate.

random.generator(1000) 

Is this approach correct?

like image 211
Lydie Avatar asked Feb 02 '16 08:02

Lydie


1 Answers

The first step is to take the inverse of your cdf function, which in this case can be done with simple arithmetic:

invcdf <- function(y) a2 * log(a1/(1-y) - 1) + a3

Now you want to call the inverse cdf with standard uniformly distributed random variables to sample:

set.seed(144)
a1 <- 1 ; a2 <- 2 ; a3 <- 3
invcdf(runif(10))
#  [1] -2.913663  4.761196  4.955712  3.007925  1.472119  4.138772 -3.568288
#  [8]  4.973643 -1.949684  6.061130

This is a histogram of 10000 simulated values:

hist(invcdf(runif(10000)))

enter image description here

And here is the plot of the pdf:

x <- seq(-20, 20, by=.01)
plot(x, df(x))

enter image description here

like image 123
josliber Avatar answered Sep 30 '22 02:09

josliber