7

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?

josliber
  • 43,891
  • 12
  • 98
  • 133
Lydie
  • 117
  • 1
  • 10

1 Answers1

4

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

josliber
  • 43,891
  • 12
  • 98
  • 133
  • thanks for your detailed answer it is very helpful. I see that you made changes both in my cdf function and the inverse that you proposed! Is there any error in my CDF expression? – Lydie Feb 02 '16 at 11:02
  • to explain more, I have initially the CDF function defined by : cdf=function(x) { 1 - a1/(1+exp((x-a3)/a2)) The density function is given by the derivate : df=function(x) { ((-a1/a2)*exp((x-a3)/a2))/(1+exp((x-a3)/a2))^2 } – Lydie Feb 02 '16 at 11:08
  • @N.Fk yes, initially your cdf took values close to 1 for low inputs and values close to 0 for high inputs, so I think you had accidentally flipped it -- I simply did one minus the expression you had before. – josliber Feb 02 '16 at 13:59
  • If I want to make some conditions in generated numbers, do you have any idea how doing this ? I mean that I want for example the values generated to be grater or equal to 100? Is that possible ? I already try to repeat the process until having values satsifying my constraints but it is not efficient! – Lydie Feb 04 '16 at 13:46
  • @Is.Fk Since that sounds like a different question, could I suggest that you ask it separately using the "Ask Question" button? That way the whole community can work to answer it. – josliber Feb 04 '16 at 15:24