Monte Carlo method and its accuracy

Metd Monte Carlo means numerical solution.
math problems using random variable simulation. An idea of the history of the method and the simplest examples of its use can be found in Wikipedia .

In the method itself there is nothing difficult. It is this simplicity that explains the popularity of this method.

The method has two main features. The first is the simple structure of the computational algorithm. The second is a calculation error, usually proportional
$\ sqrt {D \ zeta / N}$ where $D \ zeta$ - some constant, and $N$ - the number of tests. It is clear that to achieve high accuracy on this path is impossible. Therefore, it is usually said that the Monte Carlo method is especially effective in solving those problems in which the result is needed with a little accuracy.
')
However, the same problem can be solved by different versions of the Monte Carlo method, which correspond to different values $D \ zeta$ . In many problems, it is possible to significantly increase the accuracy by choosing a method of calculation, which corresponds to a significantly lower value. $D \ zeta$ .

General scheme of the method

Suppose we need to calculate some unknown quantity m. Let's try to invent such a random variable. $\ xi$ so that $M \ xi = m$ . Let at the same time $D \ xi = {{b} ^ {2}}$ .
Will consider $N$ independent random variables ${{\ xi} ^ {1}}, {{\ xi} ^ {2}}, \ ldots, {{\ xi} ^ {N}}$ (implementations) whose distributions coincide with the distribution $\ xi$ . If a $N$ is large enough, then according to the central limit theorem the distribution of the sum ${{\ rho} _ {N}} = \ sum \ limits_ {i} {{{\ xi} _ {i}}}$ will be approximately normal with parameters $M \ rho_N = Nm$ , $D \ rho_N = Nb ^ 2$ .

On the basis of the Central Limit Theorem (or, if you wish, the Muavre-Laplace limit theorem), it is not difficult to obtain the relation:

$P \ left (\ left | \ frac {{{\ rho} _ {N}}} {N} -m \ right | \ le k \ frac {b} {\ sqrt {N}} \ right) = P \ left (\ left | \ frac {1} {N} \ sum \ limits_ {i} {{{\ xi} _ {i}}} - m \ right | \ le k \ frac {b} {\ sqrt {N }} \ right) \ to 2 \ Phi (k) -1,$

Where $\ Phi (x)$ - distribution function of the standard normal distribution.

This is an extremely important relation for the Monte Carlo method. It gives the calculation method $m$ , and an error estimate.

In fact, we find $N$ random values $\ xi$ . From this ratio, it is clear that the arithmetic average of these values will be approximately equal to $m$ . With probability close to $(2 \ Phi (k) -1)$ the error of such an approximation does not exceed the magnitude $kb / \ sqrt {N}$ . Obviously, this error tends to zero with increasing $N$ .

Depending on the goals, the last relation is used in different ways:

If we take k = 3, we get the so-called “rule $3 \ sigma$ ":

$P \ left (\ left | \ frac {{{\ rho} _ {N}}} {N} -m \ right | \ le 3 \ frac {b} {\ sqrt {N}} \ right) \ approx 0.9973 .$
If a specific level of computation reliability is required $\ alpha$ ,
$P \ left (\ left | \ frac {{{\ rho} _ {N}}} {N} -m \ right | \ le \ left ({{\ Phi} ^ {- 1}} \ left ((1 + \ alpha) / 2 \ right) \ right) \ frac {b} {\ sqrt {N}} \ right) \ approx \ alpha$

Calculation accuracy

As can be seen from the above ratios, the accuracy of the calculations depends on the parameter $N$ and values $b$ - standard deviation of a random variable $\ xi$ .

At this point, I would like to indicate the importance of the second parameter. $b$ . This is best shown by example. Consider the calculation of a definite integral.

The calculation of a definite integral is equivalent to the calculation of areas, which gives an intuitively clear algorithm for calculating the integral (see the Wikipedia article). I will consider a more effective method (a special case of the formula for which, however, is also in the article from Wikipedia). However, not everyone knows that instead of a uniformly distributed random variable, this method can use almost any random variable specified on the same interval.

So, it is required to calculate a certain integral:

$I = \ int \ limits_ {a} ^ {b} {g (x) dx}$

Choose an arbitrary random variable $\ xi$ with distribution density ${{p} _ {\ xi}} (x)$ determined on the interval $(a, b)$ . And consider the random variable $\ zeta = g (\ xi) / {{p} _ {\ xi}} (\ xi)$ .

The expected value of the last random variable is:

$M \ zeta = \ int \ limits_ {a} ^ {b} {[g (x) / {{p} _ {\ xi}} (x)] {{p} _ {\ xi}} (x) dx = I}$

Thus, we get:

$P \ left (\ left | \ frac {1} {N} \ sum \ limits_ {i} {{{\ zeta} _ {i}}} - I \ right | \ le 3 \ sqrt {\ frac {D \ zeta} {N}} \ right) \ approx 0.9973.$

The last ratio means that if you choose $N$ values ${{\ xi} ^ {1}}, {{\ xi} ^ {2}}, \ ldots, {{\ xi} ^ {N}}$ then at a sufficiently large $N$ :

$\ frac {1} {N} \ sum \ limits_ {i} {\ frac {g ({{\ xi} _ {i}})} {{{p} _ {\ xi}} ({{\ xi} _ {i}})} \ approx I}$

.

Thus, to calculate the integral, you can use almost any random variable. $\ xi$ . But dispersion $D \ zeta$ , and with it the accuracy estimate, depends on which random variable $\ xi$ take for calculations.

It can be shown that $D \ zeta$ will have a minimum value when ${{p} _ {\ xi}} (x)$ proportional to | g (x) |. Choose such a value ${{p} _ {\ xi}} (x)$ In the general case, it is very difficult (complexity is equivalent to the complexity of the problem being solved), but it is worth being guided by this consideration, i.e. choose a probability distribution in the form similar to the module of the integrable function.

Numerical example

The theory, of course, is a good thing, but let's look at a numerical example: $a = 0$ ; $b = \ pi / 2$ ; $g (x) = cos (x)$ .

We calculate the value of the integral using two different random variables.

In the first case, we will use a uniformly distributed random variable on [a, b], i.e. ${{p} _ {\ xi}} (x) = 2 / \ pi$ .

In the second case, we take a random variable with a linear density on [a, b], i.e. ${{p} _ {\ xi}} (x) = \ frac {4} {\ pi} (1-2x / \ pi)$ .

Here is a graph of the specified functions.

It is easy to see that the linear density better matches the function. $g (x)$ .

The program code of the model example in the mathematical package Maple

restart; with(Statistics): with(plots): #  g:=x->cos(x): a:=0: b:=Pi/2: N:=10000: #  p1:=x->piecewise(x>=a and x<b,1/(ba)): p2:=x->piecewise(x>=a and x<b,4/Pi-8*x/Pi^2): # plot([g(x),p1(x),p2(x)],x=a..b, legend=[g,p1,p2]); #   I_ab:=int(g(x),x=0..b); #  -      #       INT:=proc(g,p,N) local xi; xi:=Sample(RandomVariable(Distribution(PDF = p)),N); evalf(add(g(xi[i])/p(xi[i]),i=1..N)/N); end proc: #   I_p1:=INT(g,p1,N);#c   I_p2:=INT(g,p2,N);#c   #  Delta1:=abs(I_p1-I_ab);#c   Delta2:=abs(I_p2-I_ab);#c   #    delta1:=Delta1/I_ab*100;#c   delta2:=Delta2/I_ab*100;#c   #  Dzeta1:=evalf(int(g(x)^2/p1(x),x=a..b)-1); Dzeta2:=evalf(int(g(x)^2/p2(x),x=a..b)-1); #     3*sqrt(Dzeta1)/sqrt(N); #     3*sqrt(Dzeta2)/sqrt(N);

The file with this program can be found here.

The exact value of the integral is easy to calculate analytically, it is equal to 1.

The results of one simulation with $N = 10$ :

For a uniformly distributed random variable: $I \ approx 1.21666$ .

For a random variable with a linear distribution density: $I \ approx 0.97641$ .

In the first case, the relative error is more than 21%, and in the second 2.35%. Accuracy $3 \ sqrt {\ frac {D \ zeta} {N}}$ in the first case it is equal 0.459, and in the second - 0.123.

I think this model example shows the importance of choosing a random variable in the Monte Carlo method. By choosing the correct random value, you can get a higher accuracy of calculations, with fewer iterations.

Of course, one-dimensional integrals are not calculated this way, there are more accurate quadrature formulas for this. But the situation changes during the transition to multidimensional integrals, since quadrature formulas become cumbersome and complex, and the Monte Carlo method is applied only with minor modifications.

Number of iterations and random number generators

It is not difficult to see that the accuracy of the calculations depends on the number $N$ random variables included in the amount. Moreover, to increase the accuracy of calculations 10 times you need to increase $N$ 100 times.

When solving some problems, it is necessary to take a very large number to obtain an acceptable accuracy of the estimate. $N$ . And given that the method often works very quickly, then implementing the latter with modern computing capabilities is not at all difficult. And the temptation is to simply increase the number $N$ .

If a physical phenomenon is used as a source of randomness (a physical sensor of random numbers), then everything works fine.

Often, pseudo-random number sensors are used for Monte-Carlo calculations. The main feature of such generators is the presence of a certain period.

The Monte Carlo method can be used with values $N$ not exceeding (preferably a lot smaller) the period of your pseudo-random number generator. The latter fact follows from the condition of independence of random variables used in modeling.

For large calculations, you need to make sure that the properties of the random number generator allow you to do these calculations. In standard random number generators (in most programming languages), the period most often does not exceed 2 in the degree of bitness of the operating system, or even less. When using such generators, one must be extremely careful. It is better to study the recommendations of D. Knut, and build your own generator, which has a known and sufficiently long period in advance.

Literature

Popular lectures on mathematics 1968. Issue 46. Sobol I.M. Monte Carlo method. M .: Science, 1968. - 64 p.

Source: https://habr.com/ru/post/274975/

All Articles