# How to decide convergence of the ensemble average in a Monte Carlo process

I am simulating electromagnetic scattering off a rough surface. The usual process is to do a Monte Carlo simulation, which is briefly described as follows.

• Generate a randomly rough surface, and compute the scattered far-field intensity in a particular direction of interest.
• Repeat this process for several different instances of rough surfaces, and get an ensemble average of the far-field intensity.

This ensemble average converges after sufficient instances have been taken. So, my question is that how does one determine in a rigourous manner if convergence has happened? In the electromagnetic scattering community, this question is not addressed, and heuristic estimates are used (e.g. "we used 100 instances") instead.

Thanks apriori.

-
In case you are not aware, there is a scientific computation site now in beta which may have a higher density of experts. I'm sure they would be happy to have this question (and it might be a better fit there than here), but (1) I won't migrate it unless you ask and (2) I don't think it is in any danger of being closed here. Your choice. –  dmckee Jul 27 '12 at 0:14
Thanks for the tip @dmckee, I will keep that group in mind for future questions. I've got one good answer here already, and I don't mind if you migrate it to the other group. –  udax Jul 27 '12 at 17:39

In this problem, you have a real valued quantity which has an unknown distribution, and you want to estimate the error on the average from averaging N independent draws. This is an application of the central limit theorem.

If you choose $N$ independent picks from a distribution $\rho(x)$ with zero mean and finite second moment, meaning that

$$\int x^2 \rho(x) = \sigma^2$$

Then add the picks together to get $X$, you get that the sum is distributed as the $N$-fold convolution of $\rho$ with itself. This eventually converges to a Gaussian of second moment $N\sigma^2$. This means that the error in taking an average over $N$ steps is the width of the distribution of $X/N$, which is

$${\sqrt{N\sigma^2}\over N} = {\sigma\over \sqrt{N}}$$

In other words, the error is Gaussian distributed with a width that falls off as the square root of the number of trials. Your error is the $\sigma$ of the distribution divided by the square root of the number of trials.

If your trials give you $x_1,x_2,...,x_N$ for x, you calculate the mean of these, and call it x:

$$x = {\sum_i x_i \over N}$$

Then calculate the average square of $(x_i-x)$ to get a best-estimator for $\sigma$

$$\sigma^2 = {\sum_i (x_i - x)^2 \over N}$$

Then your error is $\sigma/\sqrt{N}$. This is the best estimate to make in your circumstances.

The one thing you have to check is the convergence to a Gaussian. This means that you make sure that the roughness model gives a distribution of intensities where the $sigma$ is not dominated by rare events, so that the estimate from N trials of $\sigma$ is reliable. In principle, you could have a contrived rough surface model which produces ridiculous things on very rare occasions--- for example, suppose the roughness occasionally conspires over long distances to be a nearly exactly periodic sinusoid. This makes the surface a diffraction grating, with extremely sharp peaks at certain directions. If you are looking at the direction where the diffraction peak occurs, you get a huge outlier in the scattered intensity for a diffraction grating configuration, and this might occur in too few trials, so that the peak might never show up in your limited number of trials.

In practice, you just make sure your roughness model is not conspiratorial (so as to make periodic gratings) and for this it is usually enough (aside from contrived nonlocal roughness conspiracy) to check that the distribution of intensities doesn't have a power-law tail. You cannot do anything rigorous without a description of the roughness model, so that you can show that diffraction grating conspiracies don't happen with any reasonable probability, so that the second-moment is well described by what you see.

This is usually true, absent obvious tail phenomenon (for example, you see 2 trials out of 200 that dominate the variance) so doing a central limit theorem error analysis is good enough.

-
Thanks for the useful answer, @Ron. I do know the statistics of the rough surface, and can therefore be confident that there are no pathological surfaces in my ensemble. Just to be precise, what you mean by error is the standard deviation. PS: A small typo in your answer: the second moment of the sum is N \sigma^2, not N \sigma. –  udax Jul 27 '12 at 20:35

Do you have access to the standard deviation? If so, the coefficient of variation is a good measure of convergence: std.dev/mean. A generalisation is the root mean square deviation. Those wikipedia pages also tell you how to handle, to some extent, situations where the mean is close to zero.

-
Well, since I have access to the intensity for each instance of a rough surface, I can generate all statistics, including the standard deviation. std.dev/mean will fail for a zero-mean process, though. –  udax Jul 27 '12 at 17:41
@UdX I just updated the answer with some links that may be useful. Hope that helps. –  Gabriel Landi Jul 27 '12 at 23:13