5
$\begingroup$

Implement an estimator using Monte Carlo integration of the quantity $$\theta=\int_0^1e^{-x^2}(1-x)dx$$ Estimate $\theta$ with a variance lower than $10^{-4}$ by writing the variance of this estimator depending on sample size.

We can write $$\theta=\int \phi(x)f(x)dx$$ where $\phi(x)$ is a function and $f(x)$ is a density so that $$\phi(x)f(x)=e^{-x^2}(1-x)\mathbb{I}_{(0,1)}(x)$$ The exercise leaves open the choice of the density. Thus the estimator has the form $$\hat{\theta}=\frac{1}{n}\sum_i \phi(x_i)$$ The exercise asks for an estimate of $\theta$ with variance lower than $0.0001$ by expressing the variance of the estimator as a function of n.

$\endgroup$
5
  • $\begingroup$ What's wrong with just increasing $n$ to 10^7? That seems to do it. $\endgroup$ May 7, 2015 at 1:59
  • $\begingroup$ @philchalmers But how can I write the variance of my estimator in terms of n? $\endgroup$
    – user72621
    May 7, 2015 at 2:03
  • 1
    $\begingroup$ You have the variance estimator given $n$ already in your answer. MC integration error is root-n consistent, so the variance decreases to 0 at a rate of 1/$n$. $\endgroup$ May 7, 2015 at 2:53
  • 1
    $\begingroup$ you seem to have a typo for the line X = (exp(1)^(-U^2)*(1-U)): it'll be equal to (1 - U) $\endgroup$ May 7, 2015 at 9:40
  • 1
    $\begingroup$ This appears to be just another version of your question at stats.stackexchange.com/questions/151224/…, which also asks how to perform a MC integration with a given estimation variance. $\endgroup$
    – whuber
    May 8, 2015 at 13:27

3 Answers 3

6
$\begingroup$

The problem is that without knowing exactly what $\theta$ is, we cannot know the variance of its Monte-Carlo estimator. The solution is to estimate that variance and hope the estimate is sufficiently close to the truth.


The very simplest form of Monte-Carlo estimation surrounds the graph of the integrand, $f(x) = e^{-x^2}(1-x)$, by a box (or other congenial figure that is easy to work with) of area $A$ and places $n$ independent uniformly random points in the box. The proportion of points lying under the graph, times the area $A$, estimates the area $\theta$ under the graph. As usual, let's write this estimator of $\theta$ as $\hat\theta$. For examples, see the figure at the end of this post.

Because the chance of any point lying under the graph is $p = \theta / A$, the count $X$ of points lying under the graph has a Binomial$(n, p)$ distribution. This has an expected value of $np$ and a variance of $np(1-p)$. The variance of the estimate therefore is

$$\text{Var}(\hat \theta) = \text{Var}\left(\frac{AX}{n}\right) = \left(\frac{A}{n}\right)^2\text{Var}(X) = \left(\frac{A}{n}\right)^2 n \left(\frac{\theta}{A}\right)\left(1 - \frac{\theta}{A}\right) = \frac{\theta(A-\theta)}{n}.$$

Because we do no know $\theta$, we first use a small $n$ to obtain an initial estimate and plug that into this variance formula. (A good educated guess about $\theta$ will serve well to start, too. For instance, the graph (see below) suggests $\theta$ is not far from $1/2$, so you could start by substituting that for $\hat\theta$.) This is the estimated variance,

$$\widehat{\text{Var}}(\hat\theta) = \frac{\hat\theta(A-\hat\theta)}{n}.$$

Using this initial estimate $\hat\theta$, find an $n$ for which $\widehat{\text{Var}}(\hat\theta) \le 0.0001 = T$. The smallest possible such $n$ is easily found, with a little algebraic manipulation of the preceding formula, to be

$$\hat n = \bigg\lceil\frac{\hat\theta(A - \hat\theta)}{T}\bigg\rceil.$$

Iterating this procedure eventually produces a sample size that will at least approximately meet the variance target. As a practical matter, at each step $\hat n$ should be made sufficiently greater than the previous estimate of $n$ so that eventually a large enough $n$ is guaranteed to be found for which $\widehat{\text{Var}}(\hat\theta)$ is sufficiently small. For instance, if $\hat n$ is less than twice the preceding estimate, use twice the preceding estimate instead.


In the example in the question, because $f$ ranges from $1$ down to $0$ as $x$ goes from $0$ to $1$, we may surround its graph by a box of height $1$ and width $1$, whence $A=1$.

One calculation beginning at $n=10$ first estimated the variance as $2/125$, resulting in a guess $\hat n = 1600$. Using $1600$ new points (I didn't even bother to recycle the original $10$ points) resulted in an updated estimated variance of $0.0001545$, which was still too large. It suggested using $\hat n = 2473$ points. The calculation terminated there with $\hat\theta = 0.4262$ and $\widehat{\text{Var}}(\hat\theta) = 0.00009889$, just less than the target of $0.0001$. The figure shows the random points used at each of these three stages, from left to right, superimposed on plots of the box and the graph of $f$.

Figure

Since the true value is $\theta = 0.430764\ldots$, the true variance with $n=2473$ is $\theta(1-\theta)/n = 0.00009915\ldots$. (Another way to express this is to observe that $n=2453$ is the smallest number for which the true variance is less than $0.0001$, so that using the estimated variance in place of the true variance has cost us an extra $20$ sample points.)

In general, when the area under the graph $\theta$ is a sizable fraction of the box area $A$, the estimated variance will not change much when $\theta$ changes, so it's usually the case that the estimated variance is accurate. When $\theta/A$ is small, a better (more efficient) form of Monte-Carlo estimation is advisable.

$\endgroup$
6
  • $\begingroup$ That you did is "Hit or Miss"? $\endgroup$
    – user72621
    May 10, 2015 at 13:50
  • $\begingroup$ So you're … estimating the variance of the estimate? So the variance has a variance? $\endgroup$
    – endolith
    May 25, 2022 at 0:17
  • 1
    $\begingroup$ @endolith That is correct. $\endgroup$
    – whuber
    May 25, 2022 at 12:37
  • $\begingroup$ "Let's write this as $\hat\theta$." Does $\hat\theta$ = "The proportion of points lying under the graph, times the area A" or just "The [number] of points lying under the graph"? "this" is ambiguous, but I think you mean the latter $\endgroup$
    – endolith
    May 25, 2022 at 20:57
  • 1
    $\begingroup$ @endolith Thank you--that is indeed ambiguous within the context of the immediately preceding sentence. I'll fix it. The convention of putting a hat on an estimator of a quantity, though, rescues the notation from that ambiguity. $\endgroup$
    – whuber
    May 25, 2022 at 21:01
7
$\begingroup$

Implement an estimator using Monte Carlo integration of $$\theta=\int\limits_0^1e^{-x^2}(1-x)dx$$

While you can use a $\mathcal{U}([0,1])$ distribution for your Monte Carlo experiment, the fact that both $$x \longrightarrow \exp\{-x^2\}\quad \text{and}\quad x \longrightarrow (1-x)$$ are decreasing functions suggest that a decreasing density would work better. For instance, a truncated Normal $\mathcal{N}^1_0(0,.5)$ distribution could be used: \begin{align*}\theta&=\int\limits_0^1e^{-x^2}(1-x)\,\text{d}x\\&=[\Phi(\sqrt{2})-\Phi(0)]\sqrt{2\pi\frac{1}{2}}\int\limits_0^1\frac{1}{\Phi(\sqrt{2})-\Phi(0)}\dfrac{e^{-x^2/2\frac{1}{2}}}{\sqrt{2\pi\frac{1}{2}}}(1-x)\,\text{d}x\\&=[\Phi(\sqrt{2})-\Phi(0)]\sqrt{\pi}\int\limits_0^1\frac{1}{\Phi(\sqrt{2})-\Phi(0)}\dfrac{e^{-x^2}}{\sqrt{\pi}}(1-x)\,\text{d}x\end{align*} which leads to the implementation

n=1e8
U=runif(n)
#inverse cdf simulation
X=qnorm(U*pnorm(sqrt(2))+(1-U)*pnorm(0))/sqrt(2)
X=(pnorm(sqrt(2))-pnorm(0))*sqrt(pi)*(1-X)
mean(X)
sqrt(var(X)/n)

with the result

>     mean(X)
[1] 0.4307648
>     sqrt(var(X)/n)
[1] 2.039857e-05

fairly close to the true value

> integrate(function(x) exp(-x^2)*(1-x),0,1)
0.4307639 with absolute error < 4.8e-15

Another representation of the same integral is to use instead the distribution with density$$f(x)=2(1-x)\mathbb{I}{[0,1]}(x)$$and cdf $F(x)=1-(1-x)^2$ over $[0,1]$. The associated estimation is derived as follows:

> x=exp(-sqrt(runif(n))^2)/2
> mean(x)
[1] 0.4307693
> sqrt(var(x)/n)
[1] 7.369741e-06

which does better than the truncated normal simulation.

$\endgroup$
9
  • $\begingroup$ You have chosen the normal distribution, but the distribution support has to be compatible with the limits of integration? Or the support may be different since it encompasses the integration range?To write the variance in function of n, I can simply do $\sigma^2/n$ where $\sigma^2$ is the value of variance? $\endgroup$
    – user72621
    May 7, 2015 at 11:47
  • 2
    $\begingroup$ I used a truncated Normal for this purpose. Hence the normalising constant with $\Phi(1)-\Phi(0)$. $\endgroup$
    – Xi'an
    May 7, 2015 at 12:27
  • 1
    $\begingroup$ Why the variance is $\frac{\sigma}{\sqrt(n)}$ and not $\frac{\sigma^2}{n}$? $\endgroup$
    – user72621
    May 7, 2015 at 16:49
  • 1
    $\begingroup$ en.wikipedia.org/wiki/Standard_error#Standard_error_of_the_mean $\endgroup$
    – Zen
    May 8, 2015 at 13:11
  • 1
    $\begingroup$ @whuber: (+1) for your answer which addresses the second part of the question. I used 1e8 only to follow suit from earlier answers. $\endgroup$
    – Xi'an
    May 8, 2015 at 20:20
0
$\begingroup$

It's not clear whether "write the variance of estimator" means to write the equation or the results of the execution. If the latter is the case then all you need to do is to run your code at different $n$ and show how the variance shrinks with $n$.

If the former is the case, then you have to show the equation for the variance estimate of the Monte Carlo algorithm.

$\endgroup$
2
  • $\begingroup$ In reality it is to do both. $\endgroup$
    – user72621
    May 7, 2015 at 13:01
  • $\begingroup$ In this case show the estimate of the variance using MC convergence theorems. Those are usually based on some form of central limit theorem application. $\endgroup$
    – Aksakal
    May 7, 2015 at 13:14

Your Answer

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