# Importance Sampling to evaluate integral in R

I have asked the question here also. However, there might be something wrong with my theoretical understanding hence I'm asking here as it is more relevant. Kindly do not diss without looking first.

I have referred to these links - here and here. The latter link is a question of mine, but the answers, as they turn out, aren't quite right.

• I have the following integral $$\int_{0.01}^1 x^{-0.5}\,\text{d}x$$ To solve this using Importance Sampling MC integration, one needs to select an importance pdf that is approximately the same as the function plot

My R code to solve the same is this :

#function 1 - importance sampling
w <- function(x) dunif(x,0.01,1)/dbeta(x,0.7,1)
f <- function(x) x^(-0.5)
X <- rbeta(1000,0.7,1)
Y <- w(X)*f(X)
c(mean(Y),var(Y))


True integral value - 1.8

Using the Importance Sampling code above - 1.82 (where my importance PDF is Beta(0.7,1)

which is quite alright so I'm assuming the code is correct.

• Now I have this integral $$\int_{0.3}^8 [1+\sinh(2x)\log(x)]^{-1}\text{d}x$$

for which my code is :

#function 2
w <- function(x) dunif(x,0.01,1)/dnorm(x,0.5,0.25)
f <- function(x) (1+sinh(2*x)*log(x))^(-1)
X <- rnorm(1000,0.5,0.25)
Y <- w(X)*f(X)
Y <- Y[!is.na(Y)]
c(mean(Y),var(Y))


True Integral Value ~0.7014

Value from executing above code ~3.25 (where my importance PDF is normal(0.5,sd=0.25)

What am I doing wrong?

1) Take function to be evaluated as f(x).

2) Generate samples from Importance PDF g(x) that is truncated between the intervals.

3) Get the mean(f(x)/g(x)) which is the integral.

• have you looked at the plot formed by the mean for each iteration (meaning for each new random value) ? You could see if you have converged yet. 1000 random values seems pretty low to me. – Scratch Mar 6 '14 at 9:57
• @Scratch - 1000 is not a good number yes but it shouldn't be significantly different, right? For example, if 0.7 is supposed to be the answer, 0.5-0.7 is fine. I have looked at more numbers though yes and looked at the histogram. – Raaj Mar 6 '14 at 23:03
• Well if you have tried with 100 000 and the difference is still significant, the explanation is somewhere else... – Scratch Mar 7 '14 at 13:42
• No I'm saying that the difference isn't significant when I increase the sample size . So the implementation for importance sampling might be wrong but I can't put a finger on it . – Raaj Mar 7 '14 at 14:37

When running the code provided for the second function I get

> c(mean(Y),var(Y))
[1] 3.2981238 0.5203621
> integrate(f,0.01,1)
3.19264 with absolute error < 1.1e-06


which means that the true value of the integral is close to 3.2, not to 0.70.

If you want to integrate f from 0.3 to 8, then the importance function must be

> w <- function(x) dunif(x,0.3,8)/dnorm(x,0.5,0.25)


and I changed the handling of the NA's in the problem by

> f <- function(x) (x>0)*(1+sinh(2*x)*log(abs(x)))^(-1)


> Y <- w(X)*f(X)
> integrate(f,0.3,8)
2.77512 with absolute error < 2.4e-05
> integrate(f,0.3,8)$val/(8-.3) [1] 0.3604052  which shows a good fit even though the Normal importance function N(0.5,0.5) may be too concentrated, i.e., it has too small a variance to cover (0.3,8) with some reasonable probability. If one changes the scale of the Normal, > X <- rnorm(1e5,0.5,2.25) > Y <- w(X)*f(X) > c(mean(Y),var(Y)) [1] 0.3664046 1.0230197 > integrate(f,0.3,8)$va/(8-.3)
[1] 0.3604052


which shows a fairly good fit. When using a smaller standard deviation like $\sigma=0.1$, I find a poorer fit

> c(mean(Y),var(Y))
[1]   0.3370849 484.6246147


even when compared with $\sigma=1$:

> c(mean(Y),var(Y))
[1] 0.3606815 0.3232931


To make sense of those variations, I ran the experiment for a range of values of $\sigma$, from $0.1$ to $4$, with $10^6$ normal simulations, and got the following graph (with a log scale on the first axis)

which shows the improvement brought by large enough $\sigma$'s. And an optimum around $\sigma=0.74$.

Here is the relevant R code

w <- function(x) dunif(x,0.3,8)/dnorm(x,0.5,sigma)
sigs=varz=meanz=seq(.1,4,le=50)
for (i in 1:50){
sigma=sigs[i];X=0.5+sigma*X0;Y=w(X)*f(X)
varz[i]=var(Y);meanz[i]=mean(Y)}