Monte-Carlo simulation

Question

We're asked to find a monte-carlo estimate for : $$ \theta\:=\: \int^1_0 (e^x -1)dx$$

No further info is provided. And I have no idea how to proceed.

What I know is , if we are required to estimate say $E(Y)= \mu$ , we simulate $n$ values from the distribution of $Y$ , say , $ y_1 ,y_2,y_3,....,y_n$ and hence estimate $\mu$ as $\hat{\mu} = \dfrac{\sum_{i=1}^n y}{n}$.

Can anyone help me with this ?

Answer 1

Normally, you can always regard an integral as an expectation. For example, in your case, the variable $x$ ranges from $0$ to $1$ so it's reasonable to assume a beta distribution whose special case is the uniform distribution. You may have heard of the law of unconscious statistician and we're gonna use it here. Define $f\left(x\right)=e^{x}-1$. Then $$ \mathrm{E}\left(\dfrac{f\left(X\right)}{\mathcal{B}\left(X\,|\,\alpha,\beta\right)}\right)=\int_{0}^{1} \dfrac{f\left(x\right)}{\mathcal{B}\left(x\,|\,\alpha,\beta\right)} \mathcal{B}\left(x\,|\,\alpha,\beta\right)\, dx $$ which is exactly the same as your integral. By the way, $\mathcal{B}\left(x\,|\,\alpha,\beta\right)$ denotes the density function of the beta distribution with parameters $\alpha,\beta$ evaluated at $x$. So by the strong law of large numbers, generate a massive number of beta variables and compute the following. Let's assume you have generated $M$ variables. $$ \dfrac{1}{M}\sum_{m=1}^{M} \dfrac{f\left(X^{(m)}\right)}{\mathcal{B}\left(X^{(m)}\,|\,\alpha,\beta\right)} \overset{M\to\infty}{\rightarrow} \int_{0}^{1} f\left(x\right)\, dx $$ where $X^{(m)}\sim\mathcal{B}\left(\alpha,\beta\right)$ independently. However, the reason why we use uniform distribution, $\mathcal{U}\left(0,1\right)\overset{d}{\equiv}\mathcal{B}\left(1,1\right)$, is because its PDF is conveniently reduced to a single constant $1$ so the calculation gets much simpler. $$ \dfrac{1}{M}\sum_{m=1}^{M}f\left(U^{(m)}\right)\overset{M\to\infty}{\rightarrow} \int_{0}^{1}f\left(x\right)\, dx $$ where $U^{(m)}\sim\mathcal{U}\left(0,1\right)$ independently.

To summarize, you should look at the lower bound and the upper bound of the integral and pick an appropriate distribution that has the same support. Make the function into a form of an expectation which is essentially the same as the original integrand. Generate a large enough number of random variables and rely on the law of large numbers.

It will be clear, once you try both, that the results will be the same no matter which one of the above two you pick to generate random variables from. Additionally, Monte-Carlo integration converges quite slowly. So you will have to sample a lot of random numbers.

Answer 2

A good starting point is https://en.wikipedia.org/wiki/Monte_Carlo_integration.

The goal is to find the volume (surface) under the curve $$f(x)=e^x-1$$ in the interval $\Omega=[0,1]$, which is equal to $$e^1-1-e^0 \approx 0.7182.$$

Monte Carlo integration samples $n$ times from the parameter space $\Omega$ and evaluates $f_i(x_i)$ for $i=1,...n$. Subsequently the MC-integral (estimate) is given by

$$n^{-1} \sum_n f_i$$

Note that for multivariate parameter space the integral is given by

$$V n^{-1} \sum_n f_i$$

where $V = \int_\Omega dx$ grows with the number of parameters.

For your example, a simple R code solves the problem:

f = function(x) exp(x) - 1 #define function

n = 10^6               # Number of MC draws
u<-runif(n,0,1)        # Draw n times from uniform distribution in (0,1)
mean(f(u))             # Estimate the integral as simple mean
[1] 0.7184804          # Varies across draws

exp(1)-1 - exp(0)      # True value
[1] 0.7182818