CODEWITHSUNDEEP
ralgorithmsimulation
user9802913
user9802913

Reputation: 245

Simulation of random variables using acceptance-rejection method

I have the following algorithm

Step 1. Simulate the value of Y with qj=P(Y=j)

Step 2. Generate an uniform variable

Step 3. If U<=Pj/(c*qj) then X=j and stop. Otherwise come back to step 1.

And the specific example:

X=1,2,3,4 with P1=.2,P2=.15,P3=.25,P4=.4

Generate values for X

  1. Let Y~UD(1,4)

  2. c=.4/.25

Here is my approach to implement this algorithm in R:

f<-function(){
n<-4

probY<-c(.25,.25,.25,.25)
probX<-c(2,.15,.25,.4)

X<-rep(0,4)

U<-runif(n)

c<-.4/.25

for(j in 1:4)
{
if(U[j]<= probX[j]/(c*probY[j]))
   X<-j

}


return(X)

}

The output is 4 I don't think is correct. I am not sure If I should have write Y<-runif(n,1,4) nor this line probY<-c(.25,.25,.25,.25). Also this line 'Otherwise come back to step 1.' is missing in the loop, though is always the same .25

Could someone help please?

Upvotes: 1

Views: 1035

Answers (1)

BurlyPotatoMan
BurlyPotatoMan

Reputation: 194

I think the problem here is a bit of confusion with how the algorithm works.

In order to generate a single value from your distribution (X = 1,2,3,4 where P(X = 1) = .2, P(X = 2) = .15, P(X = 3) = .25, P(X = 4) = .4), we need to follow the steps of the algorithm. Assuming we choose c = .4/.25:

1.Generate y from Y ~ UD(1,4).

2.Generate u from U ~ U(0,1).

3.Check to see if u≤f(y)/cg(y). If it is, define x = y and you're done! If not, go back to step 1.

In the code you gave, you never actually generate a y variable. Here is a function that should work for you instead! Hopefully my comments explain it well enough!

accRej <- function(){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(1)
  i <- 1

  #Now, start the loop!
  while(i <= 1){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i (making it 2) to stop the while loop
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}

Notice in this code that probX[i] is equal to f(y) in our algorithm, and since Y~UD(1,4), .25 = g(y) always. Hope this helps!

Also, here is code to generate n random variables with this method. It is essentially the same as above, just with an option to change 1 to n.

accRej <- function(n){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(n)
  i <- 1

  #Now, start the loop!
  while(i <= n){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}

Upvotes: 2

Related Questions