# How to use control variate method to estimate $\theta = \int_{1}^{\infty}\frac{x^2}{\sqrt{2\pi}}e^{-x^2/2}dx$

This is the target I want to integrate with Monte Carlo with control variate method: $$\theta = \int_{1}^{\infty}\frac{x^2}{\sqrt{2\pi}}e^{-x^2/2}dx$$ I have checked with Wolfram that it is 0.400626, so the control variate should be converge to this value. I use standard normal distribution and gamma distribution (shape = 3 and rate = 1) as two control variate, but fail to converge it! What's wrong with my code? or Idea?

Here is my R code and output.

sample_size <- seq(from = 100, to = 10^4, by = 10)
target <- function(x){
x^2 * exp(-(x^2)/2) /sqrt(2*pi)
}
Sim4.1.theta <- numeric(length(sample_size))
Sim4.1.se <- numeric(length(sample_size))
Sim4.2.theta <- numeric(length(sample_size))
Sim4.2.se <- numeric(length(sample_size))
MC.sim.4 <- function(size){
u1 <- rnorm(size)
f2.1 <- u1 <- u1[u1>=1]
T1.1 <- target(u1)

u2 <- rgamma(size, shape = 3, rate = 1)
f2.2 <- u2 <- u2[u2>=1]
T1.2 <- target(u2)

c.star.1 <- -lm(T1.1~f2.1)$$coeff c.star.2 <- -lm(T1.2~f2.2)$$coeff
T2.1 <- T1.1 + c.star.1*(f2.1 - pnorm(1,lower.tail = FALSE))
T2.2 <- T1.2 + c.star.2*(f2.2 - pgamma(1,shape = 3, rate = 1, lower.tail = FALSE))

control1.estimate <- mean(T2.1[u1>=1])
control2.estimate <- mean(T2.2[u2>=1])

control1.se <- sd(T2.1)/sqrt(size)
control2.se <- sd(T2.2)/sqrt(size)
return(rbind(control1.estimate,control2.estimate,control1.se,control2.se))
}
for (i in 1:length(sample_size)) {
tem <- MC.sim.4(sample_size)
Sim4.1.theta[i] <- tem
Sim4.1.se[i] <- tem
Sim4.2.theta[i] <- tem
Sim4.2.se[i] <- tem
}
plot(x = sample_size, y = Sim4.1.theta, type = 'l',col = '#2166AC', ylim = c(0,0.5), xlab = '# of sampling size')
lines(x = sample_size, y = Sim4.2.theta, col = '#B2182B')
abline(a=0.400626,b=0,col='red') • "the control variate should be converge to this value": this is not right, a control variate is an integrand with known expectation, not equal to the unknown integral you are trying to find. Apr 11 at 6:57

Your basic Monte Carlo integral is wrong, before you even add in the control variates. When I return T1.1 from MC.sim.4 I get

> MC.sim.4(100000)
[,1]
simple.MC         0.2546318389
control1.estimate 0.3854293298
control2.estimate 0.2284484298
control1.se       0.0001006455
control2.se       0.0002170077


As you say, the true value of the target is 0.40026

> integrate(function(x) x*x*dnorm(x),lower=1,upper=Inf)
0.400626 with absolute error < 5.7e-07


But that's not what your T1.1 is doing.

> z<-rnorm(1e5)
> mean((z*z*dnorm(z))[z>1])
 0.2549178


There are two problems. First, you don't want the normal density in target, because you get the Normal density by sampling from a Normal. Second, you don't want to drop the values with z<1; you want to set them to zero

Fixing these problems

> z<-rnorm(1e5)
> mean((z*z)*(z>1))
 0.4006146


The basic code for the control variates looks ok, but you'll have to fix similar problems in how those variables are defined.