# Monte Carlo integration aim for maximum variance

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.

-

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

-
the mortal flaw in your argument in substituting the unconditional probability with conditional. When you sample a point from V, it only makes a sense to talk about $\pi_V=Pr(v\in T|v\in V)$. The probability which you are talking about is $\pi=Pr(v\in T)$, i.e. unconditional. Of course, we want $\pi$ equal to 100%, i.e. we want the enclosing volume to be as close to the volume of interest, that's trivial. However, within this enclosing volume we want the sampling points to be picked at $\Pi_V$ equal to 50%. Think about it. – Aksakal May 2 '14 at 18:01
@Aksakal You seem to be describing importance sampling. – whuber May 2 '14 at 18:43
yes. And you are saying that it's best to know the answer in advance, so there's no need to sample at all. – Aksakal May 2 '14 at 19:06
@Aksakal Yes! That's correct. Knowing the ideal to attain, though, provides useful guidance concerning how to choose the shape of the enclosing region $V$. The variance formula allows one to contemplate what can be gained in accuracy (or reduced computational effort) in return for the additional complexity of using more complicated regions $V$. For instance, I could have decreased the $.077$ variance to $.004$ by choosing the linear function $7x+2y-14$ for the top of $V$ rather than a constant function--but then I would need to find an efficient way to sample uniformly from this $V$. – whuber May 2 '14 at 19:26
Going further, it is easy to sample uniformly from the region between $7x+2y-14$ and $7x+2y-15$, which sandwiches the graph of $xy$ in the region of interest. Doing so would lower the variance to $0.000015$, implying that six-figure accuracy could be achieved with a few thousand samples rather than ten million samples (required by the smallest possible height of $21$). If I interpreted your answer correctly, this is the kind of thing you were contemplating. (With this particular solution $75\%$ of samples are expected to lie within $V$.) – whuber May 2 '14 at 19:42

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$, $$\textrm{Var}\left(\hat{p}\right) = \frac{p(1-p)}{N},$$

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|$: $$\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}.$$ Now, apply the fact that $p$ depends on $|T|$ and $V$: $$=\frac{V^2\,(|T|\,(V-|T|))/V}{N} = \frac{|T|V^2-|T|^2V}{N}.$$ 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$.

-
+1 In an effort to make this argument clearer, I have posted a reply that contains a similar analysis. – whuber May 2 '14 at 16:26

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.

-
I thought this was a really great intuitive answer. It would be helpful if whoever marked it down could provide some constructive criticism – M. Berk May 2 '14 at 15:03
@M.Berk Although I did not downvote this answer, I believe it to be incorrect. It does reflect a valuable insight, however: Monte-Carlo integration can be improved by removing a portion $T_0\subset T$ of known measure from the target region $T$ and using Monte-Carlo methods to estimate the measure of what remains, $T-T_0$. – whuber May 2 '14 at 16:24