Introduction

I would like to write a simple memo for Monte Carlo error estimation in this document. I will show matimatical formulation at first, after that, check the result with R.

Mathmatical formulation

To evaluate an integral numerically, we approximate it with Monte Carlo as the following: \[ S = \int dx f(x)p(x), \\ \hat{S} = \frac{1}{N} \sum_{i=1}^{N}f(x_i), x_i \sim p(x) \] where \(p(x)\) is a probabilistic density function satisfing the condition \(\int dx p(x) = 1\).

At first, we define \[ \bar{S} = \mathbb{E}\left[\hat{S}\right] = \frac{1}{N}\mathbb{E}\left[\sum_{i=1}^N f(x_i)\right] = \frac{1}{N}\sum_{i=1}^N \mathbb{E}\left[f(x_i)\right] = \frac{1}{N}\sum_{i=1}^N \bar{f(x)} = \bar{f(x)}. \]

To estimate Monte Carlo error \(\sigma_S\) included in \(\hat{S}\), we just need to evaluate the following equation: \[ \sigma_S^2 = \mathbb{E} \left[ \left( \hat{S} - \bar{S} \right)^2 \right] \\ = \mathbb{E}\left[\hat{S}^2\right] - \bar{S}^2 \\ = \frac{1}{N^2} \mathbb{E}\left[ \sum_{i,j=1}^{N} f(x_i)f(x_j)\right] - \bar{f(x)}^2 \\ = \frac{1}{N^2} \mathbb{E}\left[ \sum_{i=1}^{N} f(x_i)^2\right] + \frac{1}{N^2} \mathbb{E}\left[ \sum_{i,j=1, i \neq j}^{N} f(x_i)f(x_j)\right] - \bar{f(x)}^2 \\ = \frac{1}{N^2} N\mathbb{E}\left[ f(x_i)^2\right] + \frac{1}{N^2} \sum_{i,j=1, i \neq j}^{N} \mathbb{E}\left[ f(x_i)\right] \mathbb{E}\left[ f(x_j)\right] - \bar{f(x)}^2 \\ = \frac{1}{N} \mathbb{E}\left[ f(x_i)^2\right] + \frac{1}{N^2} N(N-1) \bar{f(x)}^2 - \bar{f(x)}^2 \\ = \frac{1}{N} \left( \mathbb{E}\left[ f(x_i)^2\right] - \bar{f(x)}^2 \right) \\ = \frac{1}{N} \sigma_f^2 . \] It means that we just need to evaluate sample error \(\sigma_f\) if we want to know Monte Carlo error \(\sigma_S\).

Example

If we evaluate the value \(\int_0^{L} dx x^2, L \in \mathbb{R}\), we re-write this as the following: \[ S = \int_0^{L} dx x^2 = L \int_0^{L} dx x^2 \frac{1}{L}. \] In this case, \(\frac{1}{L}\) corresponds to the probabilistic density function(uniform distribution). We can calculate the exact solution of this; \[ \left[ \frac{1}{3} x^3 \right]_0^L = \frac{1}{3}L^3. \]

we can use R to check this analysis.

L <- 3
#Monte Carlo size("N" in this document)
size <- 10^3
#Get the uniformly distributed radom number(range:[0, L])
x <- runif(size, max=L)
#Exact solution
1/3*L^3
## [1] 9
#Monte Carlo value
f <- L*x^2
mean(f)
## [1] 9.336064
#Monte Carlo error
sd(f)/sqrt(size)
## [1] 0.2522806