Monte Carlo error propagation

Question

I want to assign an error to the standard deviation computed with a Monte Carlo error propagation method.

Now, I explain better.

If we have a random variable $x$, with mean value $x_0$ and standard deviation $\Delta x$, and a function $f(x)$, we know the mean value of $f(x)$ to be $f(x_0)$, while the standard deviation of $f(x)$ can be computed through the first-order formula:

$${\displaystyle \sigma_{f}={\sqrt {\left({\left.\frac {\partial f}{\partial x}\right|_{x=x_0}} \right)^{2}\Delta x_{}^{2} }}} $$

But this formula is not really useful when $\Delta x / x_0$ is large or when $f$ is not linear around $x_0$. Hence, we can estimate the standard deviation of $f(x)$ with a Monte Carlo simulation: we generate $N$ random numbers with mean value $x_0$ and standard deviation $\Delta x$, and we apply the function $f(x)$ to each of them. Then we compute the standard deviation of $f(x)$ this way:

$$ {\displaystyle \sigma_f={\sqrt {\frac {\sum _{i=1}^{N}(f(x_{i})-{f(x_0)})^{2}}{N}}}} $$

Now $\sigma_f$ changes from a Monte Carlo simulation to another, and $\sigma_f$ is more accurate for large $N$.

My question is:

How can I evaluate the error on $\sigma_f$?

I expect this error, which is the standard deviation of $\sigma_f$, to be dependent on $N$ and $x_0/\Delta x$, and on the function $f$.

Answer 1

The easiest way to do this is to monitor the convergence rate of your standard deviation, $\sigma_f$. In general, you would run for some $N$ samples, then generate another sample $N+1$ and see how much the prediction of $\sigma_f$ changes.

In practice, to avoid noise in the convergence history, you would typically check in between batches of runs. For example, you run $N$ samples and compute $\sigma_f(N)$, then run another $M$ samples (like, say, 1000 or whatever is appropriate) and then compute $\sigma_f(N+M)$ again. And you then cutoff your simulation when the change in predicted values with a new batch is less than some threshold or flatlines.

Monte Carlo methods tend to converge rather poorly, $O(\sqrt{N})$ in general. But, they converge at that rate regardless of the number of dimensions of the uncertain space, so it avoids the curse of dimensionality that plagues other methods.