Monte Carlo integration with imposed variance

Question

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.

Answer 1

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$.

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.

Answer 2

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.

Answer 3

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.