# Estimate integral value using Monte Carlo Importance Sampling method

I have to estimate the value of this integral:

$\int_{0}^{0.5713107589} e^{-3.9365491x}dx$

using Monte Carlo Importance Sampling method.

If I understood the method correctly, to estimate the value of the integral I have to get the mean of $\frac{f(x_i)}{w(x_i)}$, with $x_i$ being generated randomly according to the weight function $w(x)$. Is this correct?

If it is, how do I figure out the function $w(x)$? And after getting the function, how can I generate random numbers according to this function?

edit: My attempt: I tried to generate a function that has the same behaviour of the one I want to integrate.

The red curve is the function I want to estimate the integral. The blue one is a beta distribution with $a = 1$ and $b = 3$.

I don't know if this is correct and I don't know what I should try now. Is my weight function correct? How can I generate random numbers from it?

edit 2:

My code in R is:

monte_carlo_importance_sampling <- function(f, n) {
xi <- rbeta(n, 1, 3)
w <- function(x) dbeta(x, 1, 3)
integral <- mean(f(xi)/w(xi))
print(integral)
}


f being the function I want to integrate.

This code gives me the correct answer for the integral from 0 to 1, and I guess it's not a coincidence, right? It's because the function rbeta(n, 1, 3) generates numbers in the interval (0, 1), am I correct? So, do I have to generate random numbers from 0 to 0.5713107589? Do you have a tip to generate such numbers?

edit 3: I followed the steps of this link: http://www.lce.hut.fi/teaching/S-114.1100/lect_9.pdf

Instead of using the beta distribution I started my work again with the function $w(x) = Ce^{-x}$

Then I used Wolfram to do $\int_0^{0.5713108} Ce^{-x}dx = 1$ and I got $C = 2.29771$ so that $w(x) = 2.29771e^{-x}$.

The next step was doing $y = \int_0^x w(x)dx = 2.29771 - 2.29771e^{-x}$ and getting the invese $x = -log(-0.435216(y-2.29771))$.

My code:

# Functions
w = function(x) 2.29771*exp(-x)
x = function(y) -log(-0.435216*(y-2.29771))
f <- function(x) exp(-3.9365491*x)

monte_carlo_importance_sampling <- function(f, n, w, x) {
xi <- x(runif(n, 0, 1))
integral <- mean(f(xi)/w(xi))
print(integral)
}


And I got the result $0.227224$, which is really close to the correct value from Wolfram, $0.227228$.

What do you think? Is my work correct? Thank you.

• Welcome to our site! This looks like a homework or textbook problem - please add the [self-study] tag & read its wiki. – Silverfish Apr 17 '16 at 22:32
• @Silverfish, thank you. I tried to improve my post. – f.c.r Apr 17 '16 at 23:31
• If I'm understanding correctly, wouldn't $e^{-3.9365491x}$ itself be the weighting function? In other words, you're trying to give preference to values of $x$ where $e^{-3.9365491x}$ is larger and thus contributes more to the total integral. Or am I misunderstanding? – barrycarter Apr 18 '16 at 3:19
• For the edit 3 solution, you may as well use the exponential with the right number, i.e., $C\exp\{−2.29771x\}$ in which case your importance sampling estimate has zero variance. – Xi'an Apr 20 '16 at 9:23

Indeed, importance sampling is about approximating the value of an integral$$\mathfrak{I}=\int_\mathcal{A} f(x)\text{d}x$$using a manageable density $\omega$ with [at least] the same support as $\mathcal{A} \cap \text{supp}(f)$ and easy to simulate. A sample simulated from $\omega$ is then produced, $x_1,\ldots,x_n$, and the average$$\frac{1}{n}\sum_{i=1}^n \dfrac{f(x_i)}{\omega(x_i)}\mathbb{I}_\mathcal{A}(x_i)$$is an unbiased estimator of $\mathfrak{I}$ with variance going down to zero as $\text{O}(n^{-1})$ [if the variance is finite].
Given that you picked a Beta B(1,3) as your choice of $\omega$, you simply need to find a random generator for this distribution. If working with R, you should check the details of the function rbeta.
• The integral can be seen as $$\int_0^1 \tilde{f}(x)\text{d}x$$where$${\tilde{f}} (x) = \mathbb{I}_{(0,0.5713107589)}(x)\exp\{ -2.29771 x\}$$ i.e., you can move the integral bounds inside the function. – Xi'an Apr 20 '16 at 9:17