Monte Carlo method is a very powerfull technique for estimating integral but is not efficient since its convergence rate is low. It typically has a variance of the form \(\frac { { \sigma }^{ 2 } }{ n }\). One method to reduce the variance is to increase the sample size (n) but this would be computationally costly. In order to achieve better estimation, variance reduction techniques were developped.
These techniques can be placed into 2 groups;
* Improve sampling methods
+ antithetic sampling
+ stratification
+ Importance sampling + common random numbers
We will outline the antithetic sampling method.
Let us consider 2 indentically distributed random variables \({ U }_{ 1 }\) and \({ U }_{ 2 }\).
In general,
\(Var(\frac { { U }_{ 1 }+{ U }_{ 2 } }{ 2 } )\quad =\quad \frac { 1 }{ 4 } (Var({ U }_{ 1 })\quad +\quad Var({ U }_{ 2 })\quad +\quad 2Cov({ U }_{ 1 },{ U }_{ 2 }))\quad\)
If \({ U }_{ 1 }\) and \({ U }_{ 2 }\) are indepent, then;
\(Var(\frac { { U }_{ 1 }+{ U }_{ 2 } }{ 2 } )\quad =\quad \frac { 1 }{ 4 } (Var({ U }_{ 1 })\quad +\quad Var({ U }_{ 2 }))\quad\)
so the variance of \(({ { U }_{ 1 }+{ U }_{ 2 }) }/{ 2 }\) is smaller if \({ U }_{ 1 }\) and \({ U }_{ 2 }\) are negatively correlated than when the variables are independent.
In antithetic sampling method, we consider U and 1-U as our uniform random distribution, this preclude that g(x) is monotone on interval [0, 1]
Let us illustrate with examples;
first, let us consider:
\(\theta =\int _{ 0 }^{ 1 }{ { e }^{ -x }dx }\)
\({ e }^{ -x }\) is monotonic on interval [0,1] so we can apply the Anithetic Sampling method.
and
\(\theta =\int _{ 0 }^{ 1 }{ \frac { 1 }{ 1+x } } dx\)
again, function \(\frac { 1 }{ 1+x }\) is monotonic on interval [0,1]
Calculations will be done in R.
f_antithetic <- function (g, n, antithetic = TRUE){
### function to compute regular Monte carlo or antithetic monte carlo
u <- runif(n/2, 0, 1)
if(!antithetic){
v <- runif(n/2, 0, 1) # u and v are independent
} else {
v <- 1-u
}
u <- c(u,v) # combine u,v
Y <- g(u)
theta_hat <- mean(Y)
return(theta_hat)
}
f_expo <- function (x){
y<-exp(-x)
return(y)
}
f_oneoverx <- function (x){
y<-1/(1+x)
return(y)
}
# for validating the results, we will also calculate the integral
f_int_exp <- function(x){
y <- -exp(-x)
return(y)
}
f_int_oneover <- function(x){
y <- log(1+x)
return(y)
}
#######
# We will repeat each calculation 4 times
multiplier <- c(1, 10, 100, 500)
r <- length(multiplier)
r_mc <- numeric()
r_anti <- numeric()
r_true <- numeric()
r_m <- numeric()
for (i in 1:r) {
# we will use based on multiplier
m <- 100*multiplier[i]
v_mc <-f_antithetic(f_expo, m, FALSE)
v_anti <-f_antithetic(f_expo, m, TRUE)
v_true <- f_int_exp(1)-f_int_exp(0)
r_m <- c(r_m, m)
r_mc <- c(r_mc, v_mc)
r_anti <- c(r_anti, v_anti)
r_true <- c(r_true, v_true)
}
print(round(rbind(r_m, r_true, r_mc, r_anti), 5))## [,1] [,2] [,3] [,4]
## r_m 100.00000 1000.00000 1.0000e+04 5.0000e+04
## r_true 0.63212 0.63212 6.3212e-01 6.3212e-01
## r_mc 0.68362 0.63628 6.3172e-01 6.3302e-01
## r_anti 0.63408 0.63216 6.3216e-01 6.3177e-01r2_mc <- numeric()
r2_anti <- numeric()
r2_true <- numeric()
r2_m <- numeric()
for (i in 1:r) {
# we will use based on multiplier
m <- 100*multiplier[i]
v_mc <-f_antithetic(f_oneoverx, m, FALSE)
v_anti <-f_antithetic(f_oneoverx, m, TRUE)
v_true <- f_int_oneover(1)-f_int_oneover(0)
r2_m <- c(r2_m, m)
r2_mc <- c(r2_mc, v_mc)
r2_anti <- c(r2_anti, v_anti)
r2_true <- c(r2_true, v_true)
}
print(round(rbind(r2_m, r2_true, r2_mc, r2_anti), 5))## [,1] [,2] [,3] [,4]
## r2_m 100.00000 1000.00000 1.0000e+04 5.0000e+04
## r2_true 0.69315 0.69315 6.9315e-01 6.9315e-01
## r2_mc 0.72109 0.69419 6.9408e-01 6.9297e-01
## r2_anti 0.69327 0.69328 6.9355e-01 6.9349e-01From the result, it is clear that we get better estimation with the antithetic sample method that the regular Monte Carlo estimation.