How to use self-normalized importance sampling method to estimate $\int_{1}^{\infty} \frac{x^2}{\sqrt{2\pi}}e^{\frac{-x^2}{2}}dx$?

Question

I want to use self-normalized importance sampling methods to estimate $$\int_{1}^{\infty} \frac{x^2}{\sqrt{2\pi}}e^{\frac{-x^2}{2}} \,dx$$ I choose exponential distribution with rate $\lambda=1$ as my importance function which is $$f(x)=e^{-x}$$ The true value of the integration is about $0.400$ but I get 0.799. The following is my R code. I follow the algorithm in page 32. http://people.sabanciuniv.edu/sinanyildirim/Lecture_notes.pdf I still can't find the error in my code.

N=10000
f=function(x){
  return((x^2)*(x>=1))
}
p=function(x){
  return(dnorm(x))
}
#((1/sqrt(2*pi))*exp((-x^2)/2))
q = function(x) {
  return(exp(-x))
}
x = rexp(N, rate =1)
theta.hat2=sum((p(x)/q(x))*f(x))/sum((p(x)/q(x)))
theta.hat2

---

Update:
I use t distribution as the important function because it has same support with the target density and I get the desired value. How can I compute the variance of this estimator?

N=10000
f=function(x){
  return((x^2))
}
p=function(x){
  return(exp((-x^2)/2))
}
q = function(x) {
  return(dt(x,df=3))
}
x = rt(N,df=3)
w_u=p(x)/q(x)
w=w_u/sum(w_u)
theta.hat2=sum(w*f(x)*(x>=1))
theta.hat2
# 0.4064571

Answer 1

A first remark is that, even when the normalising constants are both available (for p(⋅) and q(⋅), self normalised importance sampling, while biased, may produce a smaller mean squared error. For instance, take the extreme example of the integrand being constant: self normalised importance sampling then has a mean squared error of zero while regular importance sampling does not. We discuss the case for using self-normalised importance sampling in our book.

In the current situation, however, self-normalised importance sampling does not work. The reason being that the support of the importance function (i.e., exponential) is not equal to the support of the target density (i.e., normal) and hence that the expectation of the ratio of normalised densities is not equal to one:

>  x = rexp(N,1) #exponential sample
>  mean(exp(dnorm(x,l=T)+x)/sqrt(2*pi)) 
[1] 0.1995025

A second remark is that a basic rule in Monte Carlo integration is to never ever simulate zeros, i.e., to always use a proposal $q(\cdot)$ that shares the same support as the integrand $f(\cdot)p(\cdot)$. Rather than an $\mathcal E(1)$ proposal, which returns about 60% of zeroes, one should thus use a distribution supported over $(1,\infty)$, for instance an exponential $\mathcal E(1)$ drifted by one (1). The change in variance is massive:

> x=rexp(1e5,1)+1  #drifted exponential
> mean(x^2*exp(x-1-(x^2)/2)/sqrt(2*pi))
[1] 0.4008881
> var(x^2*exp(x-1-(x^2)/2)/sqrt(2*pi))
[1] 0.02471647
> x=x-1            #non-drifted exponential
> mean((x>1)*x^2*exp(x-(x^2)/2)/sqrt(2*pi))
[1] 0.4004227
> var((x>1)*x^2*exp(x-(x^2)/2)/sqrt(2*pi))
[1] 0.3433074

Answer 2

Your proposal density is exponential with unit rate, so you are generating proposal values:

$$X_1,X_2,...,X_N \sim \text{IID Exp}(1).$$

The integral in question is effectively being evaluated using the importance sampling approximation of the following form:

$$\begin{align} H &\equiv \int \limits_1^\infty \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \bigg( -\frac{x^2}{2} \bigg) \ dx \\[6pt] &= \int \limits_1^\infty \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \bigg( x - \frac{x^2}{2} \bigg) \cdot \exp( -x) \ dx \\[6pt] &= \int \limits_1^\infty \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \bigg( x - \frac{x^2}{2} \bigg) \cdot \text{Exp}( x|1) \ dx \\[6pt] &= \int \limits_0^\infty \frac{x^2}{\sqrt{2 \pi}} \cdot \exp \bigg( x - \frac{x^2}{2} \bigg) \cdot \mathbb{I}(x \geqslant 1) \cdot \text{Exp}( x|1) \ dx \\[10pt] &\approx \frac{1}{N} \sum_{i=1}^{N} \frac{x_i^2}{\sqrt{2 \pi}} \cdot \exp \bigg( x_i - \frac{x_i^2}{2} \bigg) \cdot \mathbb{I}(x_i \geqslant 1). \\[6pt] \end{align}$$

For large $N$ this approximation should give you a reasonable approximation to the integral, though you could improve the approximation by choosing a candidate distribution that is closer to the target function (in particular, one that shares the support of the integral). Below I show some R code where we approximate the integral using $N=10^5$ simulated values. This is compared with the result of numerical integration using the integrate function. The approximation result is reasonable close in this case.

#Generate proposal data
set.seed(5811387)
N <- 10^5
X <- rexp(N, rate = 1)

#Generate importance sampling approximation
INT.APPROX <- mean((X^2)*exp(X-(X^2)/2)*(X >= 1)/sqrt(2*pi))
INT.APPROX
[1] 0.4031962

#Compute integral using alternative method
f <- function(x) { (x^2)*exp(-(x^2)/2)/sqrt(2*pi) }
integrate(f, lower = 1, upper = Inf)
0.400626 with absolute error < 5.7e-07

Answer 3

Your integral is $$ \int x^2 1(x \ge 1)\frac{p(x)}{q(x)}q(x)dx $$ I'm using the notation of your code, not your $\LaTeX$.

You simulate from the exponential $q$, and evaluate $x_i^2 1(x_i \ge 1)\frac{p(x_i)}{q(x_i)}$ on each sample, and then take the mean.

mean((p(x)/q(x))*f(x))

That gives you the expected answer. Self-normalized importance sampling is not necessary here because you can evaluate the normalized version of all the densities. This code implements regular importance sampling, which is algorithm 3.1 in your reference.

---

A sufficient condition to guarantee consistency of the estimator of regular IS is $$ |f(x)|p(x) \ll q(x)|f(x)| $$ for all $x$. In your case, this is true because both target and proposal have a support that covers $\{x : x \ge 1\}$.

A sufficient condition to guarantee consistency of the estimator of self-normalized IS is

$$ p(x) \ll q(x). $$

This isn't true in your case because $q$ only covers the positive numbers, but $p$ extends across positive and negative numbers. This can be diagnosed with mean((p(x)/q(x))), which should be close to $1$ (for valid proposals), but is not.