Rejection sampling - total probability of acceptance

Question

I am given the following pdf $$f(x)=3 x^{2}, \quad 0 \leq x \leq 1$$ which i need to simulate by using rejection sampling.

I have used the following code below in R.

PDF_function <- function(x){3*x^2}

c <- PDF_function(1)

curve(PDF_function, 0, 1)
abline(a = c, b = 0, col='red')

X = runif(4500, 0, 1)
U = runif(4500, 0, 1)

pi_x <- function(x) {
  new_x = 3*x^2
  return(new_x)
}

count = 1
accept = c()

while(count <= 4500 & length(accept) < 1000){
  test_u = U[count]
  test_x = pi_x(X[count])/(c*dunif(X[count],0,1))
  if (test_u <= test_x){
    accept = rbind(accept, X[count])
    count = count + 1
  }
  count = count + 1
}

hist(accept)

mean(accept)
var(accept)

integrate(PDF_function, 0, 1)

1/(1*3)

length(accept)/count

I am asked to compare the total probability of acceptance and compare it with the theoretical one. If I understand this correctly then by computing length(accept)/count I get the estimated probability of acceptance but what is the theoretical one? I think the theoretical one is the area of the pdf function over the total area (under the envelope). If I integrate my function from 0 to 1 I get 1 (obviously) and the total area is 1 * (f(1)=3) = 3. So the theoretical is 1/3 ~= 0.333. Is this correct?

score 1 · Accepted Answer · answered Jun 04 '20 at 21:32

Yes, the theoretical value is $1/3$, which is the area under the curve over the amplified proposal distribution area, which is uniform in your case. You're double-counting the count, so if you remove the additional count increment inside the conditional statement you'll obtain empirical acceptance ratio close to $1/3$.

Rejection sampling - total probability of acceptance

1 Answers1