Editing Monte Carlo method (section)

== Simple Monte Carlo ==
Suppose one wants to know the expected value <math>\mu</math> of a population (and knows that <math>\mu</math> exists), but does not have a formula available to compute it. The simple Monte Carlo method gives an estimate for <math>\mu</math> by running <math>n</math> simulations and averaging the simulations' results. It has no restrictions on the probability distribution of the inputs to the simulations, requiring only that the inputs are randomly generated and are independent of each other and that <math>\mu</math> exists. A sufficiently large <math>n</math> will produce a value for <math>m</math> that is arbitrarily close to <math>\mu</math>; more formally, it will be the case that, for any <math>\epsilon > 0</math>, <math>|\mu - m| \leq \epsilon</math>.
<ref name=":1">{{Cite book |last=Owen |first=Art B. |url=https://artowen.su.domains/mc/Ch-intro.pdf |title=Monte Carlo Theory, Methods and Examples |date=2013 |publisher=Work in progress |publication-date=2013 |pages=15–36}}</ref>

Typically, the algorithm to obtain <math>m</math> is
 <small>''s'' = 0; 
 '''for''' ''i'' = 1 to ''n'' '''do'''
    run the simulation for the ''i''<sup>th</sup> time, giving result ''r<sub>i</sub>''</small>;
    <small>''s'' = ''s'' + ''r<sub>i</sub>'';
 '''repeat'''</small>
 <small>''m'' = ''s'' / ''n'';</small>

=== An example ===
Suppose we want to know how many times we should expect to throw three eight-sided dice for the total of the dice throws to be at least <math>T</math>. We know the expected value exists. The dice throws are randomly distributed and independent of each other. So simple Monte Carlo is applicable:

 <small>''s'' = 0; 
 '''for''' ''i'' = 1 to ''n'' '''do'''
    throw the three dice until ''T'' is met or first exceeded</small>; <small>''r<sub>i</sub>'' = the number of throws</small>;
    <small>''s'' = ''s'' + ''r<sub>i</sub>'';</small>
 <small>'''repeat'''
 ''m'' = ''s'' / ''n'';</small>

If <math>n</math> is large enough, <math>m</math> will be within <math>\epsilon</math> of <math>\mu</math> for any <math>\epsilon > 0</math>.

=== Determining a sufficiently large ''n'' ===

==== 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>

==== A formula when simulations' results are bounded ====
An alternative formula can be used in the special case where all simulation results are bounded above and below.

Choose a value for <math>\epsilon</math> that is twice the maximum allowed difference between <math>\mu</math> and <math>m</math>. Let <math>0 < \delta < 100</math> be the desired confidence level, expressed as a percentage. Let every simulation result <math>r_1, r_2, \ldots, r_i, \ldots, r_n</math> be such that <math>a \leq r_i \leq b</math> for finite <math>a</math> and <math>b</math>. To have confidence of at least <math>\delta</math> that <math>|\mu - m| < \epsilon/2</math>, use a value for <math>n</math> such that:

:<math>n\geq 2(b-a)^2\ln(2/(1-(\delta/100)))/\epsilon^2</math>

For example, if <math> \delta = 99\% </math>, then <math>n \geq 2(b - a)^2 \ln(2/0.01)/\epsilon^2 \approx 10.6(b - a)^2/\epsilon^2</math>.<ref name=":1" />