3
$\begingroup$

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.

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ Welcome to our site! This looks like a homework or textbook problem - please add the [self-study] tag & read its wiki. $\endgroup$ – Silverfish Apr 17 '16 at 22:32
  • $\begingroup$ @Silverfish, thank you. I tried to improve my post. $\endgroup$ – f.c.r Apr 17 '16 at 23:31
  • 2
    $\begingroup$ 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? $\endgroup$ – barrycarter Apr 18 '16 at 3:19
  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 20 '16 at 9:23
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I edited my orginal question with my R code. Do I have to generate random numbers from 0 to 0.5713107589? Do you have a tip to generate such numbers? $\endgroup$ – f.c.r Apr 20 '16 at 0:25
  • $\begingroup$ 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. $\endgroup$ – Xi'an Apr 20 '16 at 9:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.