Monte Carlo integration aim for maximum variance

Question

I have a question about Monte Carlo integration. As I understand it the method takes a region S of known volume V which contains the region T specified in the definite integral. $T \in S$.

Then random points in S are selected and it is checked if they belong to the volume that is the volume under the graph.

Whether a point falls in the volume or not is a binomial random variable and if 30% of the sampled points do lie in the volume under the graph then the estimated volume is 0.3V.

What I have read is that the optimize the method it is best to define S so that there is the greatest variance when sampling, i.e. when the probability of a point lying in the volume under the graph is 0.5.

This seems counter-intuitive as the greatest variation would cause the greatest sampling error and the widest confidence interval.

However when I tried to estimate the integral $\int_{x=2}^{3} \int_{y=6}^{7} xy\ dy\ dx$

I tried it in a region where the estimated volume was about 50% of the volume of S and in a region where the estimated volume was about 2% of S. I repeated the simulation several times and I was getting more accurate results with the 50% version.

I am stumped as to why this happened and why it is recommended to use a figure of around 50%. Can someone please explain why having a large variation is advantageous.

Answer 1

50% is wrong: the closer you can get to 100% the better off you are.

Let the measure of the target region $T$ be $t$ and the measure of the enclosing (or "probe") region $V$ be $v$. The chance of a uniformly random point in $V$ to lie in $T$ therefore is $t/v$. This Bernoulli distribution has variance

$$\frac{t}{v}\left(1-\frac{t}{v}\right).$$

The estimate of $t$ is based on $n$ independent uniform samples of $V$, which is therefore a Binomial$(n, t/v)$ variate with variance $n$ times greater than that of a single sample. When its outcome is $X$ the estimate will be the proportion $\hat{t}_n = v\left(\frac{X}{n}\right) = \left(\frac{v}{n}\right)X$. Its variance is

$$\text{Var}(\hat{t}_n) = \left(\frac{v}{n}\right)^2\text{Var}(X) = \left(\frac{v}{n}\right)^2 n\left(\frac{t}{v}\left(1-\frac{t}{v}\right)\right) = \frac{1}{n}t(v-t).$$

Because $V$ encloses $T$, $v\ge t$. The variance, being a linear function of $v$ in this interval, obviously is minimized at $v=t$. Ergo, the correct rule is to find an enclosing volume that is as close to $T$ as possible. Ideally, $V$ will be $T$ itself and the correct answer will be available upon taking $n=0$ samples!

---

As a check, I performed $500$ Monte-Carlo integrations of $\int_2^3\int_6^7 x y\ dy dx$ using enclosing volumes of constant heights $21$, $32.5$, and $84$ and $n=1000$ iterations per integration. For instance, here are the results of one of the integrations for $21$ (showing the target region $T$ beneath the blue surface graphing $xy$ over the square $[2,3]\times[6,7]$). The black dots fall within $T$ while the red dots, although still within $V$, fall outside $T$.

The results of these $500$ trials are

Height:   21     32.5   84
------------------------------
t/v:      77%    50%    19%
Mean:     16.258 16.235 16.280
Variance:  0.077  0.271  1.061
t(v-t)/n:  0.077  0.264  1.100

The true mean of $65/4$ was adequately estimated on average in all three situations. Their variances are comfortably close to the value of $t(v-t)/n$ previously derived.

Clearly 50% (the middle column) is not more accurate: its estimated variance of $0.271$ is almost four times worse than the estimated variance of $0.077$ achieved when the target region is $77\%$ of the probe region (left column). It would therefore take approximately $0.271/0.077 = 3.5$ times as many iterations in the $50\%$ configuration to achieve the same level of accuracy as the $77\%$ configuration. In fact, when I redid the $50\%$ calculation using $n=3500$ iterations per integral, the variance of $500$ trials was $0.065$. This does not differ significantly from $0.077$.

The import of this variance calculation is plain: lower variances mean either less computation or better accuracy (or both).

Answer 2

Partial solution that explains why 2 % is worse than 50 %, but does not arrive at the 50 % guideline.

The variance of the estimate $\hat{p}$ of the proportion $|T|/V$, \begin{equation} \textrm{Var}\left(\hat{p}\right) = \frac{p(1-p)}{N}, \end{equation}

is indeed maximized when $V=2|T|$. However, we are actually interested in the estimation error in the estimate of the volume $|T|$, not the proportion. Let us compute the variance of the estimate of $|T|$: \begin{equation} \textrm{Var}\left(\hat{|T|}\right) = \textrm{Var}\left(V \hat{p}\right) = V^2 \textrm{Var}\left(\hat{p}\right) = \frac{V^2 p(1-p)}{N}. \end{equation} Now, apply the fact that $p$ depends on $|T|$ and $V$: \begin{equation} =\frac{V^2\,(|T|\,(V-|T|))/V}{N} = \frac{|T|V^2-|T|^2V}{N}. \end{equation} Now, for constant sample size $N$ and true volume $|T|$, the variance of the estimator is a second-degree polynomial of $V$ which is increasing at $V\geq |T|$. However, the optimal solution based on this line of thought would be $V = |T|$!

Intuitively, if $S$ is huge, the proportion of samples in $T$ can be (in absolute terms) estimated pretty accurately to be about 0, but that does not tell us much about the volume of T. On the other hand, if we are able to set $S=T$ and know $V$, then we also know $|T|$, and thus this would be optimal.

I don't know where the 50 % guideline comes from, as the derivation in my answer would suggest to use the smallest possible $S$ that satisfies the conditions i) $S$ covers $T$, ii) we know the volume $|S|=V$ iii) we can sample points in $S$.

Answer 3

It's easy. Let's say T is a unit circle, and S is a square than contains it. When you sample from S, you want to pick points which are closer to where the circle's bound is, right? If you sample point in the center of a square, you know that they're going to be inside the circle too. There's very little information gained from checking whether (0,0) is inside the circle, in fact there's no information at all: we know it's inside the circle. And the same for (1,1): we know that it's definitely outside.

So, it makes a sense to go for points farther from the center, e.g. (0.7,0.6) - is this inside or outside the circle? Picking this point and checking will bring useful results. That was the intuition: you want to sample from the regions where the boundary of T goes, and in these regions the probabilities will tend to be far from both 0 and 1, and the farthest you can get is 0.5

UPDATE:

Less intuitive but precise answer is that your sampling distribution should be proportional to the integrand to minimize the variance. Which the same as saying that you want to sample more often from regions where your volume boundary goes through.