# x <- seq(-10, 10, by=0.01)
# plot(x, (sin(x)/x)^2, type="l")
# abline(h=0)

N <- 100000
sigma <- 2

Xj <- 0.1
X <- numeric(N)
for (j in 1:N) {
  Y <- Xj + rnorm(1, sd = sigma)
  U <- runif(1)

  alpha <- min(1, Xj^2*sin(Y)^2/Y^2/sin(Xj)^2)

  if (U <= alpha) {
    Xj <- Y
  }
  X[j] <- Xj
}

hist(X, prob = TRUE, breaks=100, xlim=c(-20,20))
x <- seq(-10, 10, by=0.01)
lines(x, 1/pi*(sin(x)/x)^2, col="red", lwd=2)

plot(X[1:100], type="l")