Estimate integral value using Monte Carlo Importance Sampling method

Question

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.

Answer 1

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.