Rejection Sampling

(a)

Rejection sampling for X,

We consider the proposal the absolutely value of standard normal, |X|, where \(X\) follows the standard normal distribution i.e., \[ g(x) = \frac{2}{\sqrt{2\pi}} \exp(-x^2/2). \]

Then \[ \frac{f(x)}{g(x)} =\sqrt{2\pi} x \exp(-x^2/2). \]

Note that, the function \(x\exp(-x^2/2)\) is in increasing function when \(0<x<1\) and an deceasing function when \(x\ge 1\). Thus, \[ \sqrt{2\pi} x \exp(-x^2/2)\le \sqrt{2\pi}\exp(-1/2)=1.520347=:M. \]

For Y

We consider the proposal uniform distribution U(0,1). Then \[ f(x)/g(x) =f_Y(y) = 30(y^2-2y^3+y^4)=30(y(1-y))^2\le 30*(1/2*(1-1/2))^2=30/16=1.875 \]

(b)

Sample X

f <- function(x){
  2*x*exp(-x^2)
}
M = sqrt(2*pi)*exp(-1/2)

n = 10000
sample_x <- rep(NA,n)
for ( i in 1:n){
  x = abs(rnorm(1,0,1))
  u = runif(1,0,1)
  while (u>f(x)/(M*2*dnorm(x))){
      x = abs(rnorm(1,0,1))
      u = runif(1,0,1)
  }
  sample_x[i] = x
}

hist(sample_x,probability = TRUE)
x = seq(0,3.5,0.001)
lines(x,f(x),col="red")

Sample Y

Sample X

f <- function(x){
  30*(y^2-2*y^3+y^4)
}
M = 1.875

sample_y <- rep(NA,n)
for ( i in 1:n){
  y = runif(1,0,1)
  u = runif(1,0,1)
  while (u>f(y)/(M)){
      y = runif(1,0,1)
      u = runif(1,0,1)
  }
  sample_y[i] = y
}

hist(sample_y,probability = TRUE)
y = seq(0,1,0.001)
lines(y,f(y),col="red")

Approximate the expectation

mean(pmax(sample_x,sample_y))

## [1] 0.9317092