CODEWITHSUNDEEP
rconditional-statementssampling
Lydie
Lydie

Reputation: 117

Conditioned random generating variables from a distribution function

My question is related to my previous one Generate random variables from a distribution function using inverse sampling Now I want to generate random variables from a distribution function using inverse sampling but the sampling should be conditioned. For example, if the inverse of my cdf is :

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

I used inverse sampling to generate 10 rv as follows :

invcdf(runif(10))

Now, the problem is that I want the values generated greater or less than a value. How should I introduce this condition in random generator?

When I use this to have value greater than 500 :

invcdf(runif(10,500,1e6))

I get this error message : Warning message: In log((a0/y) - 1) : NaNs produced

I already try to repeat the process until having values satsifying my constraints but it is not efficient!

 repeat{
   x=invcdf(runif(1))
     if(x>100){
     break
}

Upvotes: 1

Views: 236

Answers (2)

Severin Pappadeux
Severin Pappadeux

Reputation: 20090

As @spf614 noted, you'd better have checks in your function like

invcdf <- function(y) {
    if (a1 > y) {
        return( a2 * log(a1/y - 1) + a3 )
    }
    NaN
}

Then it works for all arguments

Sampling would be 

low <- ...
r <- invcdf(runif(low, a1, 1e6))

UPDATE

checking for NaNs in output

nof_nans <- sum(is.nan(r))
if (nof_nans > 0) {
    ....

Upvotes: 2

spf614
spf614

Reputation: 52

The reason that you're getting NaNs is that R is trying to take the logarithm of a negative number. Do you want the log term to be log((a1/y)-1) or log(a1/(y-1))? You currently have the function written the first way, and when you get a very high value for y, the term a1/y approaches zero (the speed with which it approaches zero depends on the value of a1). Thus, subtracting 1 gives you a negative number inside the log function. So if the term is meant to be how you have it written (log(a1/y-1)), you simply won't be able to calculate that above certain values of y.

The simple fix is just 

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

Upvotes: 1

Related Questions