# ST3456: Modern Statistical Methods 2
# Bernardo Nipoti
# Lab 4 - 16/10/18
# Material adapted from Jones, Maillardet, Robinson (2014)
# "Introduction to Scientific Programming and Simualation Using R"
###################################################################
### Part 1: Variance reduction via antithetic variables
## 1a. define the function g(x)=1-x^2 on [0,1].
g <- #...
## 1b. N=5000 realisations of the crude MC estimator, based on n=50
# uniforms (vectorized code)
N <- 5000
n <- 50
U_1 <- #... # matrix of n x N of uniforms
G_1 <- #... # matrix of the same size as U_1 obtained by evaluating g at the elements of U_1
theta_1 <- #... # mean of G_1 column by column ot obtain N realisations of theta_1
## 1c. Evaluate the variance of theta_1
# and plot a histogram of the N realisations of theta_1
var_1 <- #...
hist(theta_1,breaks=20,freq=F,xlab="theta",ylab="frequency",main="crude Monte Carlo")
## 1d. N=5,000 realisations of the antithetic estimator, each one based on
# n=50 uniforms (vectorized code)
U_a <- #... # matrix of n x N of uniforms
G_1a <- #... # matrix of the same size as U_a obtained by evaluating g at the elements of U_a
G_2a <- #... # matrix of the same size as U_a obtained by evaluating g at the elements of 1-U_a
theta_a <- #...
## 1e. Evaluate the variance of theta_a
# and plot a histogram of the N realizations of theta_a.
var_a <- #...
hist(theta_a,breaks=20,freq=F,xlab="theta",ylab="frequency",main="antithetic estimator")
## 1f. Compute the variance reduction. Is it independent of n and (approximately) equal to 93.75%?
reduction <- #...
cat("Variance of theta_1 is", var_1, "\n")
cat("Variance of theta_a is", var_a, "\n")
cat("Variance reduction is", format(reduction, dig=3), "percent\n")
# overlapping the two histograms (not requested in the lab instructions)
mt1 <- min(theta_1)
Mt1 <- max(theta_1)
hist(theta_a,breaks=20,freq=F,xlab="theta",ylab="frequency",main="crude vs antithetic",xlim=c(mt1,Mt1),col="red")
hist(theta_1,breaks=20,freq=F,xlab="theta",ylab="frequency",main="antithetic estimator",xlim=c(mt1,Mt1),col=rgb(0, 1, 0, 0.5),add=TRUE)
### Part 2: geometric interpretation of the antithetic estimator
## 2a. Adapt the program in 1b to int_0^1 h(x) dx via crude Monte Carlo,
# based on n uniforms, and evaluate the variance of the estimator.
h <- #... # definition of h, starting from g
N <- 5000
n <- 50
U_1 <- #... # matrix of n x N of uniforms
H_1 <- #... # matrix of the same size as U_1 obtained by evaluating h at the elements of U_1
theta_h1 <- #... # mean of H_1 column by column ot obtain N realisations of theta_1 for h
var_h1 <- #...
cat("Variance of theta_h1 is", var_h1, "\n")
# Is it the same as the variance computed at point 1e?
var_h1/var_a
## 2b. Run the codes for crude MC and antithetic estimator (point 1b and 1d)
# for g(x)=sin(pi x) and check the variance reduction (is it 0%?)
g_bis <- #...
# crude
U_1 <- #... # matrix of n x N of uniforms
G_1 <- #... # matrix of the same size as U_1 obtained by evaluating g_bis at the elements of U_1
theta_1_bis <- #... # mean of G_1 column by column ot obtain N realisations of theta_1
# antithetic
U_a <- #... # matrix of n x N of uniforms
G_1a <- #... # matrix of the same size as U_a obtained by evaluating g_bis at the elements of U_a
G_2a <- #... # matrix of the same size as U_a obtained by evaluating g_bis at the elements of 1-U_a
theta_a_bis <- #...
# variances
var_1_bis <- #...
var_a_bis <- #...
reduction_bis <- #...
cat("Variance of theta_1_bis is", var_1_bis, "\n")
cat("Variance of theta_a_bis is", var_a_bis, "\n")
cat("Variance reduction is", format(reduction_bis, dig=3), "percent\n")
## 2c. Run the codes for crude MC and antithetic estimator (point 1b and 1d)
# for g(x)=1-x and check the variance reduction (is it 100%?)
g_ter <- #...
# crude
U_1 <- #... # matrix of n x N of uniforms
G_1 <- #... # matrix of the same size as U_1 obtained by evaluating g_bis at the elements of U_1
theta_1_ter <- #... # mean of G_1 column by column ot obtain N realisations of theta_1
# antithetic
U_a <- #... # matrix of n x N of uniforms
G_1a <- #... # matrix of the same size as U_a obtained by evaluating g_bis at the elements of U_a
G_2a <- #... # matrix of the same size as U_a obtained by evaluating g_bis at the elements of 1-U_a
theta_a_ter <- #...
# variances
var_1_ter <- #...
var_a_ter <- #...
reduction_ter <- #...
cat("Variance of theta_1_bis is", var_1_ter, "\n")
cat("Variance of theta_a_bis is", var_a_ter, "\n")
cat("Variance reduction is", format(reduction_ter, dig=3), "percent\n")