Can a Monte Carlo simulation converge?

[s-soberh]

Hello,

I'm currently performing Monte-Carlo-Simulations by using Mathcad Prime 4.0. When I change the number of events, I always get different values. Should an Monte-Carlo-Simulation not converge with a higher number of events to a certain value?

[ValeryOchkov]

Can you show your Mathcad Prime 4 sheet?

[s-soberh]

Hello,

on Monday I could share my program because I'm not at the university at the weekend. But maybe it helps you that the process is exactly the same as described in the help (monte-carlo example).

[s-soberh]

Hi,

here is my sheet. The problem is that I always get different values depending on the number of events.

Explanation for the example:
I want to know the probability of failure from the functional context. This is defined by the area under the curve in the interval [-unlimited; 0]. But by varying the number of events, the values do not converge. Thus, depending on the calculation, I get different values.

[Werner_E]

Sure you get different values every time - its a montecarlo.

But I sure see some sort of convergence when you compare the results of increasing number of simulations.

When plotting the histogram data you may consider dividing the ordinate data by the number of runs for better comparability.

Why do you use the histogram data and not the original data to determine mean and standard deviation of the simulation data?

I attach two old simulations made with MC15 - they show the classical calculation of pi by choosing random points in a square and comparing the number of points in the inscribed circle with the total number.

You can see some sort of "convergence" but even with 10^7 points still distinct deviations from the value you would expect the simulation to "converge" to.

[s-soberh]

Unfortunately, I do not know how to record the histogram data. And I only have the simulated histogram data and no original ones.

[Werner_E]

---

@s-soberh wrote:

Unfortunately, I do not know how to record the histogram data. And I only have the simulated histogram data and no original ones.

---

Not sure what you mean.

You have the "original" simulated data in the variable Output and IMHO you should use this data to calculate mean and standard deviation.

You use only the first column of the output of "histogram" to calculate those values and so you are just using the midpoints of the intervals for your calculations and you ignore completely the frequency of the values. I am pretty sure that thats not what you intend.

BTW, you don't need to define Gauss yourself - you may use the built-in function "pnorm":

[s-soberh]

Can you give me an example?

[Werner_E]

Here is how I would do it.

I prefer putting the calculation into function which can later quick be called with various different arguments.

In the attached sheet (Prime 4 format) I had set it up with and without the "montecarlo" function.

The only argument of the simulation functions in the sheet is the number of runs, but you may make a, b, and the two means and standard deviations additional arguments if you like to play around with different values here and want to compare them side by side.

And here is what i sure would call some sort of convergence:

[s-soberh]

Thank you Werner, you are an expert.

One last question I have: If I have no sum in the function, but a multiplication, how do I analytically calculate the resulting standard deviation. See pictures below.

[Werner_E]

---

@s-soberh wrote:

Thank you Werner, you are an expert.

One last question I have: If I have no sum in the function, but a multiplication, how do I analytically calculate the resulting standard deviation. See pictures below.

---

Careful! The product of two independent normal distributed random variables is NOT normal distributed!

See for example here: https://math.stackexchange.com/questions/101062/is-the-product-of-two-gaussian-random-variables-also-a-gaussian

The consequence is that you can't calculate the probability for negative values with pnorm as I did in case of the sum. I guess you'll have to define the appropriate double integral for the resulting PDF (probability density function) - you'll find the necessary formula in the literature.

Nevertheless we can more easily calculate mean an standard deviation of the resulting distribution:

[s-soberh]

Hi,

maybe you can help me again. Can I also output a histogram of f1, f2 and f? Maybe you can make the program even more elegant, but it works

[Werner_E]

I am not an expert in these things but it may be that your calculations are wrong.

1) You calculated the standard deviation of the product of three normal distributed random variables by calculating the stdev of the product of two of them with the formula I provided and then using the same formula as product of this combined variable with the third. But the combined variable is NOT normal distributed and so it might not be valid to apply this formula!!

2) The same applies when you calculate sigma=sqrt(sigma1^2+sigma2^2). This formula works only for normal distributed random variables. The product of three is NOT normal distributed

3) And again when you define f and F. These formulas for PDF apply only to normal distributions but your is NOT. The integral might be much more complicated in your case and you may study some literature on the subject to get it right, if you really need those exact formulas.

In short the problem I see is that you combine various normal distributed random variables in various ways and treat the resulting random variable as if it would be normal distributed. But unless you combine them by simple addition the results are NOT normal distributed and the formulas you use don't fit.

[DJF]

It should, yes.  It can often take 10000 to get a good distribution shape for some problems.  Although I do question Mathcad's random number generator sometimes.  I know it's gotten stuck before and I had to shut it down to get fresh values.  Also playing with the histogram interval setting can change the apparent shape.

[LucMeekes]

NO! Certainly not!
Convergence is not a term that I think is used with stochastic processes. You can get convergence by iteration. With increasing the number of runs in Monte Carlo you will more likely get Divergence.
Example: if you throw dice over and over again, the same numbers from 1 to 6 will reappear, no matter how often you throw. It's not that when you've thrown dice 1 million times, the number three will stand out for example.
If you have a system with 10 parameters that each in some way determine an output value, and each parameter has a certain distribution of randomness, the more often you calculate the output value from randomly chosen parameters, each within its own distribution, the more likely you will find output values closer to the edge. So output values spread out, rather than coming closer.
What you will get, if you plot a histogram of the output values, is a better approximation of the distribution shape.
But note that the more input parameters you have and the higher the number of runs, the closer the distribution of the output values will tend to a normal=Gaussian distribution. And the perfect Gaussian distribution can reach all values from minus to plus infinity. Divergence, not convergence!
On the other hand it will allow you to appromate the expected value of the distribution more accurately, just by calculating the mean...
So if you're looking for the expected value, increasing the number of runs may be a good idea. But never expect the output values to converge to a number.


Success!
Luc

[DJF]

Maybe convergence isn't the right mathematical terminology, but as I interpret the question, I think it's a question on if most monte carlos will 'converge' to a stable distributed shape.  In that case the answer is usually yes.  With 5 runs the shape is erratic, with 1E6 it usually isn't.  Similarly, a monte carlo is famously used to solve for pi, and the more runs you perform the closer you get - so it can converge. 

And since items are bound in reality (e.g. dimensions on a print) actual divergence to infinity is something I don't need to consider when wearing my engineering hat.