Editing Monte Carlo method (section)

==== General formula ====
Let <math> \epsilon = |\mu - m| > 0 </math>. Choose the desired confidence level – the percent chance that, when the Monte Carlo algorithm completes, <math>m</math> is indeed within <math>\epsilon</math> of <math>\mu</math>. Let <math>z</math> be the <math>z</math>-score corresponding to that confidence level.

Let <math>s^2</math> be the estimated variance, sometimes called the “sample” variance; it is the variance of the results obtained from a relatively small number <math>k</math> of “sample” simulations. Choose a <math>k</math>; Driels and Shin observe that “''even for sample sizes an order of magnitude lower than the number required, the calculation of that number is quite stable.''"<ref name=":2">{{Cite journal |last1=Driels |first1=Morris R. |last2=Shin |first2=Young S. |date=April 2004 |title=Determining the number of Iterations for Monte Carlo Simulations of Weapon Effectiveness |url=https://apps.dtic.mil/sti/citations/ADA423541 |journal=Naval Postgraduate School Technical Report |issue=March 2003 - March 2004 |pages=10–11}}</ref>

The following algorithm computes <math>s^2</math> in one pass while minimizing the possibility that accumulated numerical error produces erroneous results:<ref name=":1" />
 <small>''s<sub>1</sub>'' = 0; 
 run the simulation for the first time, producing result ''r''<sub>1</sub>;
 ''m''<sub>1</sub> = ''r''<sub>1</sub></small>; <small>//''m<sub>i</sub>'' is the mean of the first ''i'' simulations</small>
 <small>'''for''' i = 2 to ''k'' '''do'''</small>
    <small>run the simulation for the ''i''<sup>th</sup> time, producing result ''r<sub>i</sub>'';</small>
    <small>''δ<sub>i</sub>'' = ''r<sub>i</sub>'' - ''m<sub>i</sub>''<sub>−1</sub>;</small>
    ''<small>m<sub>i</sub> = m<sub>i-1</sub></small>'' <small>+ (1/''i'')''δ<sub>i</sub>'';</small>
    ''<small>s<sub>i</sub> = s<sub>i-1</sub></small>'' <small>+ ((''i'' - 1)/''i'')(''δ<sub>i</sub>'')<sup>2</sup>;</small>
 <small>'''repeat'''
 ''s<sup>2</sup>'' = ''s<sub>k</sub>''/(''k'' - 1);</small>

Note that, when the algorithm completes, <math>m_k</math> is the mean of the <math>k</math> results.

The value <math>n</math> is sufficiently large when

:<math>n \geq s^2 z^2/ \epsilon^2.</math><ref name=":1" /><ref name=":2" />

If <math>n \leq k</math>, then <math>m_k = m</math>; sufficient sample simulations were done to ensure that <math>m_k</math> is within <math>\epsilon</math> of <math>\mu</math>. If <math>n > k</math>, then <math>n</math> simulations can be run “from scratch,” or, since <math>k</math> simulations have already been done, one can just run <math>n - k</math> more simulations and add their results into those from the sample simulations:

 <small>''s'' = ''m<sub>k</sub>'' * ''k''; 
 for i = ''k'' + 1 to ''n'' do</small>
    <small>run the simulation for the ''i''<sup>th</sup> time, giving result ''r<sub>i</sub>''</small>;
    <small>''s'' = ''s'' + ''r<sub>i</sub>'';
 ''m'' = ''s'' / ''n'';</small>