I have the function which I am trying to use calculate monte carlo estimate on using rejection sampling.
$f(x) = \int_2^\infty x^4 \sin(\pi(x+1)) e^{-x^2/2}dx$
My initial thought was that I could first limit the bounds for intergration to $\int_2^{10}$ since the function decreases to 0 quickly after x=4. Then I thought I could form a box around the function around the function and generate uniform points from that box and reject if they fall outside of the function. The rejection algorithm would looks something like
- generate $x\sim U(2,10)$
- Generate $y\sim U(-2,.5)$ (encloses min and max of function on the given range)
- accept if y < f(x) or else go back to 1
However I am running into issues here since the function does crosses the y axis, I'm not sure what I am doing is correct. After I am running 40000 samples, I get an mean of 6.19, when I know the true value is about -0.85.
Is the approach I am taking somewhat correct, or am I thinking about this wrong. Below is my code and a plot I made of the rejection / accept region based on the simulations that seems wrong to me since its not actually capturing the true function.
Mp = .5
Mm = -2
B=40000
x =runif(B,2,10)
y=runif(B,Mm,Mp)
acc = y <=(x^4)*(sin(pi*(x+1)))*exp(-x^(2)/2)
mean(x[acc])
var(x[acc])
plot(x, y, pch=".", col="red")
points(x[acc], y[acc], pch=".")
