Problem 5.14

Obtain a Monte Carlo estimate of \[\int_{1}^{\infty}\frac{x^2}{\sqrt{2\pi}}e^{-x^2/2}dx\] By Importance sampling

The candidates for the importance funcitons are only \(f_1\) and \(f_2\)

g <- function(x) {
    return((x^2)/(sqrt(2*pi))*exp(-(x^2/2))*(x>1))
}

#container for estimates
theta.hat <- se <- numeric(2)

m <- 10000

# f1(x) = e^-x, 0 > x > Inf
x <- rexp(m)
fg <- g(x)/exp(-x)
theta.hat[1] <- mean(fg)
se[1] <- sd(fg)

# f2(x) = (1+x^2)/pi, -Inf > x > Inf
x <- rcauchy(m)
fg <- g(x)/dcauchy(x)
theta.hat[2] <- mean(fg)
se[2] <- sd(fg)

rbind(theta.hat, se)
##                [,1]      [,2]
## theta.hat 0.4033587 0.3945016
## se        0.5883303 0.9515870

It looks like \(f_1\) gives us the lowest varience. Comparing the estimates to R’s estimate:

# Integrating fom 1 to Inf
integrate(g,1,Inf)
## 0.400626 with absolute error < 5.7e-07

Our approximations are faitly close, within ~0.005 from R’s estimate.