Discussion

Attempt to complete either 5.9, 5.13, or 5.14 in Rizzo’s Statistical Computing in R.

Problem 5.9

The Rayleigh density [156, (18.76)] is

\(f(x) = \frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}\) , \(x \geq 0\), \(\sigma > 0\)

Implement a function to generate samples from a Rayleigh() distribution, using antithetic variables. What is the percent reduction in variance of \(\frac{x+x'}{2}\) compared with \(\frac{x_1+x_2}{2}\) for independent \(X_1\), \(X_2\)?

Rizzo (pg 131)

Rayleigh <- function(x, m = 10000, antithetic=TRUE){
  u <- runif(m/2)
  if (!antithetic)
    v <- runif(m/2)
  else
    v <- 1 - u
  
  #calculate inv-cdf and varianace
  invcdf <- x * sqrt(-2 * log(v))
  invcdf
}

set.seed(123)
X1 <- Rayleigh(1.95, antithetic=FALSE)
set.seed(321)
X1prime <- Rayleigh(1.95, antithetic=FALSE)
set.seed(123)
X2 <- Rayleigh(1.95, antithetic=TRUE)

result1 <- (var(X1) - var(X2))/var(X1)
result2 <- (var(X1) - var(X1prime))/var(X1)

The variance reduction between \(X_1\) and \(X_2\) is -0.0150735 and the variance between \(X_1\) and \(X^{prime}_1\) is -0.0260725.

rev1 <- var(X1)/2 + var(X2)/2 + cov(X1, X2)
rev2 <- var(X1)/2 + var(X1prime)/2 + cov(X1, X1prime)
result3 <- (rev1 - rev2)/rev1 * 100

The percent reduction between \(\frac{x+x'}{2}\) and \(\frac{x_1+x_2}{2}\) is 0.2408354%.

Problem 5.13

Find two importance functions \(f_1\) and \(f_2\) that are supported on (1, \(\infty\)) and are ‘close’ to

\(g(x) = \frac{x^2}{\sqrt{2\pi}}e^{-x^2/2}\), \(x > 1\)

Which of your two importance functions should produce the smaller variance in estimating:

\(\int_1^\infty\frac{x^2}{\sqrt{2\pi}}e^{-x^2/2}dx\)

by importance sampling? Explain.

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.