Enhancing Monte-Carlo convergence (crude method)

Question

I am currently doing a project involving Monte-Carlo method. I wonder if there is papers dealing with a "learning" refinement method to enhance the MC-convergence, example :

Objective : estimate of $E(X)\thickapprox \sum _{i=1}^{10 000}X_i$

-> step 1 (500 simulations): $approx_1=\sum _{i=1}^{500}X_i$

(i) Defining and 'Acceptance interval'

$ A_1 = [approx_1-\epsilon_1,approx_1+\epsilon_1]$

where $\epsilon_1$ could be a function of the empirical variance and other statistics

-> step 2 (500 other simulations): "throwing" all simulation out of the interval $A_1$ , $approx_2=\sum _{i=1}^{500}X_i^{(2)}$

New 'acceptance interval'

$A_2 = [approx_2-\epsilon_2,approx_2+\epsilon_2]$

where $\epsilon_2 < \epsilon_1$

...

$ \epsilon \xrightarrow {} 0$

Thanks you for your help

Answer 1

If your variable of integration is truly one-dimensional, as you seem to be saying, then you should be using quadrature to evaluate the expectation integral. The computational efficiency of quadrature is much higher than Monte Carlo in one dimension (even accounting for modified sampling).

If your problem is actually multidimensional, your best bet is to use the first few iterations (you suggest 500 above) to help choose a scheme for importance sampling. Your windowing scheme is a different trick sometimes labeled stratified sampling, and tends to get tricky from a coding perspective.

To perform importance sampling, you will modify the distribution with what is known as an equivalent measure of your random samples so that most of them fall in the "interesting" region. The easiest technique is to ensure you are sampling from the multivariate normal distribution, and then shift the mean and variance of your samples such that, say, 90% of them fall within your "interesting" region.

Having shifted your samples, you need to then track their likelihood ratio (or Radon-Nikodym derivative) versus the original distribution, because your samples now need to be weighted by that ratio in your Monte Carlo sum. In the case of a shift in normal distributions, this is fairly easy to compute for each sample $\vec{x}$, as

$$ \frac{1}{\sqrt{2\pi\ \text{det}A^{-1}}}\exp{\left[-\frac12 (x-b)A(x-b)^t\right]} $$

where $A$ and $b$ are the covariance matrix and mean of your change to the original multivariate normal.