Question

Suppose we want to estimate S=E[X²] when X has the density that is proportional to q(x) = exp{-|x|³/3}.

1. Use rejection sampling to estimate E[X²]. Describe your algorithm first. Attach your R code and output. Be sure to count and report your acceptance ratio

Answer 1

For a random variable X with un-normalized density q(x) = exp(−|x|
3/3), consider
estimating σ
2 = E(X2
). By numerical integration, I found that the “true” value is
σ
2 = 0.776.
(a) We can do importance sampling using a standard normal proposal—the cubic
in the exponent of q means the tails are lighter than Gaussian which has a
quadratic in the exponent. That is, we sample Y ∼ g = N(0, 1), compute
importance weights w
?
(Y ) = q(Y )/g(Y ), and estimate σ
2 with
PN
i=1 w
?
(Yi)Y
2
P
i
N
i=1 w?
(Yi)
.
Based on a sample of size N = 1000, I got an estimate of 0.734 for σ
2
.
(b) We can employ the accept–reject strategy to sample from q, using the same
standard normal proposal as in (a). It’s not difficult to show (using calculus)
that
exp(−|x|
3
/3) ≤ A exp(−|x|
2
/2),
where A = exp(1/6) ≈ 1.18. Run accept–reject until a sample X1, . . . , XN of
size N is available, and estimate σ
2 with
1
N
Xn
i=1
X
2
i
.
Based on a sample of size N = 1000, I got an estimate of 0.728 for σ
2
.
(c) Using the sample X1, . . . , XN of size N = 1000 from part (b) above, I get 0.759
as an estimate of σ
2 based on the Philippe–Robert strategy.
(d) To compare the estimators in parts (b) and (c), I repeat the above experi-
ments 1000 times, each based on N = 1000, and compute the empirical stan-
dard deviation of the two estimators. In this case, I found that the standard
deviations for the accept–reject and Philippe–Robert methods are 0.029 and
0.009, respectively. Therefore, the Philippe–Robert method is significantly
more accurate than a basic accept–reject in this problem.

R code:-

rnorm.newton <- function(n, mean=0, sd=1) {
X <- numeric(n)
U <- runif(n)
for(i in 1:n) {
f <- function(x) pnorm(x) - U[i]
df <- function(x) dnorm(x)
X[i] <- newton(f, df, 0)$solution # newton function from my website
}
return(mean + sd * X)
}
rnorm.bm <- function(n=1, mean=0, sd=1) {
if(n %% 2 == 0) p <- n / 2 else p <- (n %/% 2) + 1
U <- runif(p)
V <- runif(p)
R <- sqrt(-2 * log(V))
Theta <- 2 * pi * U
X <- R * cos(Theta)
Y <- R * sin(Theta)
out <- mean + sd * c(X, Y)
return(out[1:n])
}
M <- 50000
o.nwtn <- system.time(nwtn <- rnorm.newton(M))[3]
o.bm <- system.time(bm <- rnorm.bm(M))[3]
o.rn <- system.time(rn <- rnorm(M))[3]
print(cbind(Newton=o.nwtn, BoxMuller=o.bm, rnorm=o.rn))
Problem 2.
rgamma.ar <- function(n, shape, scale=1) {
s <- shape
s.int <- floor(s)
b <- s / s.int
M <- gamma(s.int) / gamma(s) * b**s.int * (s * exp(-1))**(s - s.int)
f <- function(y) dgamma(y, shape=s)
Mg <- function(y) M * dgamma(y, shape=s.int, rate=1 / b)
acpt <- 0
total <- 0

X <- numeric(n)
while(acpt < n) {
total <- total + 1
Y <- sum(-b * log(runif(s.int)))
if(runif(1) <= f(Y) / Mg(Y)) {
acpt <- acpt + 1
X[acpt] <- Y
}
}
return(list(X=scale * X, rate.true=1 / M, rate.obs=acpt / total))
}
o <- rgamma.ar(1000, shape=5.5)
print(o[-1])
ylim <- c(0, 1.05 * dgamma(4.5, shape=5.5))
hist(o$X, freq=FALSE, col="gray", border="white", xlab="x", ylim=ylim, main="")
curve(dgamma(x, shape=5.5), add=TRUE)
Problem 3.
H.naive <- function(z, M) mean(rnorm(M) + rexp(M) <= z)
H.rb <- function(z, M) mean(pexp(z - rnorm(M)))
reps <- 5000
o.naive <- o.rb <- numeric(reps)
for(r in 1:reps) {
o.naive[r] <- H.naive(1, 5000)
o.rb[r] <- H.rb(1, 5000)
}
print(cbind(sd.naive=sd(o.naive), sd.rb=sd(o.rb)))
Problem 4.
qq <- function(x) exp(-abs(x)**3 / 3)
sigma2.is <- function(N) {
X <- rnorm(N)
w <- qq(X) / dnorm(X)
return(sum(w * X**2) / sum(w))
}
qq.ar <- function(N) {
X <- numeric(N)

r <- 1
A <- exp(1 / 6)
R <- function(x) qq(x) / A / exp(-x**2 / 2)
repeat {
Y <- rnorm(1)
if(runif(1) <= R(Y)) {
X[r] <- Y
r <- r + 1
}
if(r > N) break
}
return(X)
}
sigma2.pr <- function(N, X=NULL) {
if(is.null(X)) X <- sigma2.ar(N)
XX <- sort(X)
w <- (XX[-1] - XX[-N]) * qq(XX[-N])
return(sum(w * XX[-N]**2) / sum(w))
}
o.ar <- o.pr <- numeric(1000)
for(i in 1:1000) {
X <- qq.ar(1000)
o.ar[i] <- mean(X**2)
o.pr[i] <- sigma2.pr(1000, X)
}
print(cbind(AR=sd(o.ar), PR=sd(o.pr)))
6