My favorite method of variance reduction found in Rizzo text would be using Antithetic variables. The idea behind this method is to reduce variance by making a set of dependent outcomes; these dependent outcomes are made through negative correlating your input.

For example; if we had a uniformily distributed RV from [0,1]

\(U_1 = [.1 , .3 , .5]\)

\(U_{2} = [.9 , .7 , .5]\)

or

\(1-u\)

Without the antitheticals, we would use this formula to estimate the variance of two identically distributed random variables.

With the antitheticals(making variables dependent), we are able to use this formula to estimate variance of the two; which normally produces a significantly smaller variance. ( This is not always the case )

Lets look at Rizzo’s use case of antithetics.
The CDF of the standard normal distribution, is denoted as phi.
This function uses antithetical variables optionally; so we can see the difference between using them, and doing a plain old MC.

MC.Phi <- function(x, R = 10000, antithetic = TRUE) { 
  u = runif(R/2)
if (!antithetic){
  v = runif(R/2)} 
else
  {
    v = 1 - u
  }
  
  #When activating antitheticals; half of the vector is replaced by 1-vector
u = c(u, v)
cdf = numeric(length(x))

for(i in 1:length(x))
  {
    g      = x[i] * exp(-(u * x[i])^2 / 2) # Our integrand
    cdf[i] = mean(g) / sqrt(2 * pi) + 0.5 
   } 
cdf
}

x =  seq(.1, 2.5, length=5)

Phi = pnorm(x) 
set.seed(155)
MC = MC.Phi(x, anti = FALSE)
set.seed(155)
Anti = MC.Phi(x)

Here are the estimates

print(round(rbind(x, MC, Anti, Phi), 5))
##         [,1]    [,2]    [,3]    [,4]    [,5]
## x    0.10000 0.70000 1.30000 1.90000 2.50000
## MC   0.53983 0.75815 0.90370 0.97207 0.99432
## Anti 0.53983 0.75821 0.90382 0.97178 0.99288
## Phi  0.53983 0.75804 0.90320 0.97128 0.99379

We see the antithetic estimates become lower after ~ 2.5

To approximate variance reduction for any given x; Rizzo compares classic MC integration with antithetical variables.

m = 1000 
MC1 = MC2 = numeric(m) 
x = 1.95 
for (i in 1:m) {
  #Compute with and without antitheticals
MC1[i] = MC.Phi(x, R = 1000, anti = FALSE) 
MC2[i] = MC.Phi(x, R = 1000)
}
print(sd(MC1))
## [1] 0.006769195
print(sd(MC2))
## [1] 0.0004707065
print((var(MC1) - var(MC2))/var(MC1))
## [1] 0.9951647

Thats a 99.5% reduction in variance.

99.5% is massive; I don’t see a reason not to use this method over a crude monte carlo. I could definitely see using this approach to measure something like server utilization or AvG entities in system during a simulation.