# 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$?

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


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))
 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))
 0.4008881
> var(x^2*exp(x-1-(x^2)/2)/sqrt(2*pi))
 0.02471647
> x=x-1            #non-drifted exponential
> mean((x>1)*x^2*exp(x-(x^2)/2)/sqrt(2*pi))
 0.4004227
> var((x>1)*x^2*exp(x-(x^2)/2)/sqrt(2*pi))
 0.3433074

• I know how to use importance sampling but I have never thought that the exponential distribution can be drifted by one. Does random number rexp(1e5,1)+1 follow exponential distribution? If I really want to try self-normalized importance sampling in this case, is it possible? Thanks! Apr 6 at 11:31
• Quite obviously, rexp(1e5,1)+1 shifts the exponential variate by one, making the new rv supported by $(1,\infty)$. For self-normalised importance sampling, you have to find an alternative since dividing by the sum of the weights does not work. Apr 6 at 11:41
• If I choose t-distribution as an important function, will self-normalized importance sampling work? Both target density and important density have same support. My simulation result looks reasonable. Apr 6 at 12:17
• Yes, as I mentioned in my answer, it does not work only for different supports. Apr 6 at 12:41
• Nice explanation. In effect, then, it's not a Normal distribution that is being used: it is a half-Normal distribution, and thereby originates the factor of $2$ difference in the result.
– whuber
Apr 6 at 13:56

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
 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


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.

• I still can't understand why you say self-normalized importance sampling is not appropriate here. Can you give some examples? Thanks! Apr 4 at 17:09
• Because you can evaluate the normalized versions of all the densities, you don't need to self-normalize. Apr 4 at 17:14
• What does "the distribution is normailized" mean? Apr 5 at 3:50
• @Justin it means it integrates to one Apr 5 at 13:09
• The problem is not underflow: the supports of $p(\cdot)$ and $q(\cdot)$ differ, hence the normalising constant (or ratio thereof) is not properly estimated. Apr 6 at 7:50