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

Question

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[2]
  c.star.2 <- -lm(T1.2~f2.2)$coeff[2]
  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[1]
  Sim4.1.se[i] <- tem[2]
  Sim4.2.theta[i] <- tem[3]
  Sim4.2.se[i] <- tem[4]
}
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')

Answer 1

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])
[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))
[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.