# Monte-Carlo simulation

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 ?

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.

• Nice answer (+1), maybe you should be more specific when saying ''However, the reason why we use uniform distribution is because its PDF is conveniently reduced to a single constant 1 so the calculation gets much simpler'' because that only holds for the uniform distribution on [0,1] (which is the case for this question). And if you want to realy make it perfect you might add the formulas for the $1/M \sum_m f(X_m)$ where $x_m$ are draws from a uniform distribution as well ? – user83346 Sep 4 '16 at 14:39
• @fcop thanks for the recommendation! I will fix my answer accordingly. – Daeyoung Lim Sep 4 '16 at 14:44
• What is this law of the unconscious statistician?! – Xi'an Sep 5 '16 at 4:11
• @Xi'an it's when we unconsciously use the PDF of a random variable $X$ where we should actually derive the PDF of $Y=f(X)$ which is in reality a different random variable. – Daeyoung Lim Sep 5 '16 at 4:16

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

• Given that the question had the tag self-study, this (nice) answer may be too complete to help the OP solve the question by oneself! – Xi'an Sep 5 '16 at 4:13
• @Xi'an Yeah -- when I answered the question the self-study tag had not been assigned yet (I suggested it myself, but had already provided an answer when I had that thought). I guess that - although I give a numerical answer - the student may still be required to find a proper notation for a MC-estimate, perhaps combining mine and DaeyoungLim's answer – tomka Sep 5 '16 at 9:24