3
$\begingroup$

I am following the Physically Based Rendering: From Theory To Implementation book. In the exercises from Chapter 2 (Monte carlo integration) they ask to write a program to compute numerical integration of a function using traditional numerical integration techniques and Monte Carlo.

I chose the $f(x) = x^2$ to be integrated in the $[0,1]$ interval. For the default Monte Carlo Estimator $$ F_n = \frac{b-a}{n}\sum^N_{i=1}f(X_i) $$ I have no issue computing the integral. I sample $X_i \sim U(0,1)$ and evaluate the estimator.

However, in the book they explain that if $X_i \sim p(x)$ where $p(x)$ is different from a uniform then the estimator becomes $$ F_n = \frac{1}{n}\sum^N_{i=1}\frac{f(X_i)}{p(X_i)} $$ In my case I am assuming that $X_i \sim N(2,5)$. However I am not able to obtain the correct result. Do I need to discard the samples $N(2,5)$ not in $[0,1]$ ? This would reduce the total number of samples compared with the uniform case

$\endgroup$
1
  • 1
    $\begingroup$ Yes, you do. You are integrating over $[0,1]$, so points outside that interval don't count. You can always increase the sample size in the second case to compensate for the fact that well over half the samples will be discarded. $\endgroup$
    – jbowman
    Commented Dec 30, 2023 at 0:57

1 Answer 1

3
$\begingroup$

This method of integration is called importance sampling, and yes, you discard the samples outside the support of the integral under this method. However, rather than thinking of this as an additional step, you should see this as an application of the stated rule. Your integrand (written with the whole real line as the domain) is the function $f: \mathbb{R} \rightarrow \mathbb{R}$ given by:

$$f(x) = x^2 \cdot \mathbb{I}(0 \leqslant x \leqslant 1) \quad \quad \quad \text{for all } x \in \mathbb{R}.$$

Consequently, taking the sampling density $p(x) = \text{N}(x|2,5)$ and generating the resulting samples $X_1,...,X_N$ gives the resulting integral approximation:

$$\begin{align} F_N &= \frac{1}{N} \sum_{i=1}^N \frac{f(X_i)}{p(X_i)} \\[6pt] &= \frac{1}{N} \sum_{i=1}^N \frac{X_i^2 \cdot \mathbb{I}(0 \leqslant X_i \leqslant 1)}{\text{N}(X_i|2,5)} \\[6pt] &= \frac{\sqrt{10 \pi}}{N} \sum_{i=1}^N X_i^2 \exp \Bigg( \frac{(X_i-2)^2}{10} \Bigg) \cdot \mathbb{I}(0 \leqslant X_i \leqslant 1) \\[6pt] &= \frac{\sqrt{10 \pi}}{N} \sum_{i \in \mathscr{U}_N} X_i^2 \exp \Bigg( \frac{(X_i-2)^2}{10} \Bigg), \\[6pt] \end{align}$$

where $\mathscr{U}_N \equiv \{ i=1,...,N | 0 \leqslant X_i \leqslant 1 \}$ is the subset of sample indices for sample values in the integral range (i.e., discarding all the samples outside that range).


Computational implementation: In the code below we implement this calculation in R to confirm that the rule works. First we use the method of using all samples but setting the summation term to zero in the case where the generated value falls outside the support of the integrand. Using $N=10^7$ samples we obtain a result that is correct to three decimal places.

#Set parameters
N    <- 10^7
MEAN <- 2
VAR  <- 5

#Generate integrand
f <- function(x) { (x^2)*(x >= 0)*(x <= 1) }

#Generate integral approximation and true value
set.seed(1)
X <- rnorm(N, mean = MEAN, sd = sqrt(VAR))
TERMS  <- f(X)/dnorm(X, mean = MEAN, sd = sqrt(VAR))
INT.APPROX <- sum(TERMS)/N
INT.TRUE   <- 1/3

#Compare approximation with true integral
INT.APPROX
[1] 0.3331653
INT.TRUE
[1] 0.3333333

As an alternative, we also use the method of using only those samples where the generated value falls outside the support of the integrand. Using $N=10^7$ samples we obtain a result that is correct to three decimal places.

#Set parameters
N    <- 10^7
MEAN <- 2
VAR  <- 5

#Generate integrand (without restriction on the range)
f <- function(x) { x^2 }

#Generate integral approximation and true value
set.seed(1)
X <- rnorm(N, mean = MEAN, sd = sqrt(VAR))
USE    <- (X >= 0)*(X <= 1)
SUM    <- 0
if (sum(USE) > 0) {
for (i in which(USE > 0)) {
  SUM <- SUM + f(X[i])/dnorm(X[i], mean = MEAN, sd = sqrt(VAR)) } }
INT.APPROX <- SUM/N
INT.TRUE   <- 1/3

#Compare approximation with true integral
INT.APPROX
[1] 0.3331653
INT.TRUE
[1] 0.3333333
$\endgroup$
2
  • $\begingroup$ So basically you embed in your function the condition that samples outside the support should contribute 0 to the integrand. Am I correct? $\endgroup$
    – ElPotac
    Commented Dec 31, 2023 at 0:03
  • $\begingroup$ Yes, that is the way I have done it here, but if you were to discard those samples in a separate step that would also yield the correct answer. $\endgroup$
    – Ben
    Commented Dec 31, 2023 at 1:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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