Namespaces
Variants
Actions

Difference between revisions of "Monte-Carlo method"

From Encyclopedia of Mathematics
Jump to: navigation, search
(Importing text file)
 
m (tex encoded by computer)
 
Line 1: Line 1:
 +
<!--
 +
m0648701.png
 +
$#A+1 = 132 n = 0
 +
$#C+1 = 132 : ~/encyclopedia/old_files/data/M064/M.0604870 Monte\AAnCarlo method,
 +
Automatically converted into TeX, above some diagnostics.
 +
Please remove this comment and the {{TEX|auto}} line below,
 +
if TeX found to be correct.
 +
-->
 +
 +
{{TEX|auto}}
 +
{{TEX|done}}
 +
 
''method of statistical trials''
 
''method of statistical trials''
  
Line 4: Line 16:
  
 
==Simulation by random variables with given distributions.==
 
==Simulation by random variables with given distributions.==
As a rule such a simulation is achieved by a transformation of one or more independent values of a random number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648701.png" />, uniformly distributed in the interval <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648702.png" />. The sequence of  "sample"  values of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648703.png" /> is usually obtained on a computer by number-theoretic algorithms, of which the most widely used is the so-called method of residues, for example, in the form
+
As a rule such a simulation is achieved by a transformation of one or more independent values of a random number $  \alpha $,  
 +
uniformly distributed in the interval $  ( 0 , 1 ) $.  
 +
The sequence of  "sample"  values of $  \alpha $
 +
is usually obtained on a computer by number-theoretic algorithms, of which the most widely used is the so-called method of residues, for example, in the form
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648704.png" /></td> </tr></table>
+
$$
 +
u _ {0}  = 1 ,\ \
 +
u _ {n}  \equiv  u _ {n-} 1 5  ^ {2p+} 1 \
 +
(  \mathop{\rm mod}  2  ^ {m} ) ,\ \
 +
\alpha _ {n}  = u _ {n} \cdot 2  ^ {-} m .
 +
$$
  
Here <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648705.png" /> is the order of the mantissa of the computer and
+
Here m $
 +
is the order of the mantissa of the computer and
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648706.png" /></td> </tr></table>
+
$$
 +
= \max \{ {q } : {5  ^ {2q+} 1 < 2  ^ {m} } \}
 +
.
 +
$$
  
Numbers of this type are called pseudo-random numbers; they are used in statistical testing and in solving typical problems (see [[#References|[2]]]–[[#References|[6]]]). The length of the period of the above version of the method of residues is <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648707.png" />. Physical generators, tables of random numbers and quasi-random numbers are also used in the Monte-Carlo method. There are Monte-Carlo methods with a small number of playing parameters (see [[#References|[7]]]).
+
Numbers of this type are called pseudo-random numbers; they are used in statistical testing and in solving typical problems (see [[#References|[2]]]–[[#References|[6]]]). The length of the period of the above version of the method of residues is $  2  ^ {m-} 2 $.  
 +
Physical generators, tables of random numbers and quasi-random numbers are also used in the Monte-Carlo method. There are Monte-Carlo methods with a small number of playing parameters (see [[#References|[7]]]).
  
The standard method for simulating a discrete random variable with distribution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648708.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m0648709.png" /> is as follows: Put <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487010.png" /> if the chosen value of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487011.png" /> satisfies
+
The standard method for simulating a discrete random variable with distribution $  {\mathsf P} \{ \xi = x _ {k} \} = p _ {k} $,  
 +
$  k = 0 , 1 \dots $
 +
is as follows: Put $  \xi = x _ {m} $
 +
if the chosen value of $  \alpha $
 +
satisfies
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487012.png" /></td> </tr></table>
+
$$
 +
\sum _ { k= } 0 ^ { m- }  1
 +
p _ {k}  \leq  \alpha  < \
 +
\sum _ { k= } 0 ^ { m }
 +
p _ {k} .
 +
$$
  
The standard method for simulating a continuous random variable (sometimes called the method of inverse functions) is to use the, easily verified, representation <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487013.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487014.png" /> is the distribution function with given density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487015.png" />. Sometimes [[Randomization|randomization]] of the simulation is useful (in other words, the method of superposition), based on the expression
+
The standard method for simulating a continuous random variable (sometimes called the method of inverse functions) is to use the, easily verified, representation $  \xi = F ^ { - 1 } ( \alpha ) $,  
 +
where $  F $
 +
is the distribution function with given density $  f $.  
 +
Sometimes [[Randomization|randomization]] of the simulation is useful (in other words, the method of superposition), based on the expression
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487016.png" /></td> </tr></table>
+
$$
 +
f ( x)  = \sum _ { k } p _ {k} f _ {k} ( x) ;
 +
$$
  
here one first chooses a number <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487017.png" /> with distribution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487018.png" />, and then obtains a sample value <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487019.png" /> from the distribution with density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487020.png" />. In other methods of randomization certain parameters of a deterministic method for solving the problem are considered as random variables (see [[#References|[7]]]–[[#References|[9]]]).
+
here one first chooses a number m $
 +
with distribution $  {\mathsf P} \{ m = k \} = p _ {k} $,  
 +
and then obtains a sample value $  \xi $
 +
from the distribution with density $  f _ {m} $.  
 +
In other methods of randomization certain parameters of a deterministic method for solving the problem are considered as random variables (see [[#References|[7]]]–[[#References|[9]]]).
  
Another, more general, method for simulating a continuous random variable is the method of exclusion (method of selection), at the basis of which lies the following result: If <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487021.png" /> is uniformly distributed in a domain <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487022.png" />, then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487023.png" />. In the method of exclusion, choose a point <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487024.png" /> uniformly in a domain <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487025.png" /> and put <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487026.png" /> if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487027.png" />; otherwise repeat the selection of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487028.png" />, etc. For example, if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487029.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487030.png" />, one can take <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487031.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487032.png" />. The average number of operations in the method of exclusion is proportional to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487033.png" />.
+
Another, more general, method for simulating a continuous random variable is the method of exclusion (method of selection), at the basis of which lies the following result: If $  ( \xi , \eta ) $
 +
is uniformly distributed in a domain $  G = \{ {( x , y ) } : {0 \leq  y \leq  g ( x) } \} $,  
 +
then $  f _  \xi  ( x) = g ( x) / \overline{G}\; $.  
 +
In the method of exclusion, choose a point $  ( \xi _ {0} , \eta ) $
 +
uniformly in a domain $  G _ {1} \supset G $
 +
and put $  \xi = \xi _ {0} $
 +
if $  ( \xi _ {0} , \eta ) \in G $;  
 +
otherwise repeat the selection of $  ( \xi _ {0} , \eta ) $,  
 +
etc. For example, if $  a \leq  \xi \leq  b $
 +
and $  g ( x) = c f ( x) \leq  R $,  
 +
one can take $  \xi _ {0} = a + ( b - a ) \alpha _ {1} $,  
 +
$  \eta = R \alpha _ {2} $.  
 +
The average number of operations in the method of exclusion is proportional to $  \overline{G}\; _ {1} / \overline{G}\; $.
  
For many random variables, special representations of the form <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487034.png" /> have been obtained. For example, the random variables
+
For many random variables, special representations of the form $  \xi = \phi ( \alpha _ {1} \dots \alpha _ {n} ) $
 +
have been obtained. For example, the random variables
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487035.png" /></td> </tr></table>
+
$$
 +
\sqrt {- 2  \mathop{\rm ln}  \alpha _ {1} }
 +
\cdot \cos  2 \pi \alpha _ {2} \ \
 +
\textrm{ and } \ \
 +
\sqrt {- 2  \mathop{\rm ln}  \alpha _ {2} }
 +
\cdot \sin  2 \pi \alpha _ {2}  $$
  
have standard normal distributions and are independent; the random variable <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487036.png" /> has the gamma-distribution with parameter <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487037.png" />; the random variable <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487038.png" /> is distributed with density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487039.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487040.png" />; the random variable <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487041.png" /> has the beta-distribution with parameters <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487042.png" /> (see [[#References|[3]]]–[[#References|[6]]]).
+
have standard normal distributions and are independent; the random variable $  \mathop{\rm ln} ( \alpha _ {1} \dots \alpha _ {n} ) $
 +
has the gamma-distribution with parameter $  n $;  
 +
the random variable $  \max ( \alpha _ {1} \dots \alpha _ {n} ) $
 +
is distributed with density $  n x  ^ {n-} 1 $,  
 +
0 \leq  x \leq  1 $;  
 +
the random variable $  \mathop{\rm exp} \{ \sum _ {k=} 1  ^ {n} (  \mathop{\rm ln}  \alpha ) / ( p + k + 1 ) \} $
 +
has the beta-distribution with parameters $  p , n $(
 +
see [[#References|[3]]]–[[#References|[6]]]).
  
The standard algorithm for simulating a continuous random vector <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487043.png" /> is to successively choose the values of its components from conditional distributions corresponding to the representation
+
The standard algorithm for simulating a continuous random vector $  \\vec{xi} = ( \xi _ {1} \dots \xi _ {n} ) $
 +
is to successively choose the values of its components from conditional distributions corresponding to the representation
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487044.png" /></td> </tr></table>
+
$$
 +
f  {} _ {\\vec{xi} }  ( x _ {1} \dots x _ {n} )  = \
 +
f _ {1} ( x _ {1} ) f _ {2} ( x _ {2} \mid  x _ {1} ) \dots f _ {n} ( x _ {n} \mid  x _ {1} \dots x _ {n-} 1 ) .
 +
$$
  
The method of exclusion extends to the multi-dimensional case without change; it is only necessary, in its formulation, to regard <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487045.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487046.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487047.png" /> as vectors. A multi-dimensional normal vector can be simulated by using a special linear transformation of a vector of independent standard normal random variables. Special methods have also been developed for the approximate simulation of stationary Gaussian processes (see, for example, [[#References|[3]]], [[#References|[6]]]).
+
The method of exclusion extends to the multi-dimensional case without change; it is only necessary, in its formulation, to regard $  \xi $,  
 +
$  \xi _ {0} $
 +
and $  x $
 +
as vectors. A multi-dimensional normal vector can be simulated by using a special linear transformation of a vector of independent standard normal random variables. Special methods have also been developed for the approximate simulation of stationary Gaussian processes (see, for example, [[#References|[3]]], [[#References|[6]]]).
  
 
If, in a Monte-Carlo method calculation, random variables defined by a real phenomenon are simulated, then the calculation is a direct simulation (imitation) of this phenomenon. Computer simulations have been worked out for: the processes of transport, scattering and reproduction of particles: neutrons, gamma-quanta, photons, electrons, and others (see, for example, [[#References|[11]]]–[[#References|[18]]]); simulations of the evolution of ensembles of molecules for the solution of various problems in classical and quantum statistical physics (see, for example, [[#References|[10]]]–[[#References|[18]]]); simulation of queueing and industrial processes (see, for example [[#References|[2]]], [[#References|[6]]], [[#References|[18]]]); simulation of various random processes in technology, biology; etc. (see [[#References|[18]]]). Simulation algorithms are usually carefully developed; for example, they tabulate complicated functions, modify standard procedures, etc. None-the-less, a direct simulation often cannot provide the necessary accuracy in the estimates of the required variables. Many methods have been developed for increasing the effectiveness of a simulation.
 
If, in a Monte-Carlo method calculation, random variables defined by a real phenomenon are simulated, then the calculation is a direct simulation (imitation) of this phenomenon. Computer simulations have been worked out for: the processes of transport, scattering and reproduction of particles: neutrons, gamma-quanta, photons, electrons, and others (see, for example, [[#References|[11]]]–[[#References|[18]]]); simulations of the evolution of ensembles of molecules for the solution of various problems in classical and quantum statistical physics (see, for example, [[#References|[10]]]–[[#References|[18]]]); simulation of queueing and industrial processes (see, for example [[#References|[2]]], [[#References|[6]]], [[#References|[18]]]); simulation of various random processes in technology, biology; etc. (see [[#References|[18]]]). Simulation algorithms are usually carefully developed; for example, they tabulate complicated functions, modify standard procedures, etc. None-the-less, a direct simulation often cannot provide the necessary accuracy in the estimates of the required variables. Many methods have been developed for increasing the effectiveness of a simulation.
  
 
==Monte-Carlo algorithms for estimating multiple integrals.==
 
==Monte-Carlo algorithms for estimating multiple integrals.==
Suppose it is required to estimate an integral <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487048.png" /> with respect to the Lebesgue measure in an <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487049.png" />-dimensional Euclidean space <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487050.png" /> and let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487051.png" /> be a probability density such that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487052.png" /> can be written as:
+
Suppose it is required to estimate an integral $  J = \int h ( x)  d x $
 +
with respect to the Lebesgue measure in an $  s $-
 +
dimensional Euclidean space $  X $
 +
and let $  f _  \xi  ( x) $
 +
be a probability density such that $  J $
 +
can be written as:
 +
 
 +
$$
 +
J  =  \int\limits _ { X } f _  \xi  ( x)
 +
 
 +
\frac{h ( x) }{f _  \xi  ( x) }
 +
\
 +
d x  =  {\mathsf E} \zeta ,
 +
$$
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487053.png" /></td> </tr></table>
+
where  $  \zeta = h ( \xi ) / f _  \xi  ( \xi ) $.  
 +
By computer simulation of  $  \xi $
 +
it is possible to obtain  $  N $
 +
sample values  $  x _ {1} \dots x _ {N} $.  
 +
By the law of large numbers,
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487054.png" />. By computer simulation of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487055.png" /> it is possible to obtain <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487056.png" /> sample values <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487057.png" />. By the law of large numbers,
+
$$
 +
J  \approx  J _ {N}  = \
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487058.png" /></td> </tr></table>
+
\frac{1}{N}
 +
\sum _ { k= } 1 ^ { N }
  
Simultaneously it is possible to estimate the mean-square error in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487059.png" />, that is, the quantity <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487060.png" />, and to approximately construct a suitable confidence interval for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487061.png" />. By the choice of the density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487062.png" /> it is possible to arrange for estimates with possibly smaller variance. For example, if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487063.png" />, then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487064.png" /> and if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487065.png" />, then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487066.png" />. The corresponding algorithm is called essential sampling (choice by importance). Another common modification — the method of selection of principal part — occurs when a function <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487067.png" /> is determined whose integral is known. It is sometimes useful to combine Monte-Carlo methods and classical quadratures (cf. [[Quadrature|Quadrature]]) in so-called random quadrature formulas, the basic idea of which is that the nodes and coefficients in any quadrature sum (for example, in interpolation) are chosen randomly from a distribution, providing an [[Unbiased estimator|unbiased estimator]] of the integral [[#References|[3]]]. Particular cases of these formulas are: the so-called method of stratified sampling, in which the nodes are chosen one in each part of a fixed partition of the domain of integration and the coefficients are proportional to the corresponding volumes; and the so-called method of symmetric sampling, which, in the case of integration over <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487068.png" />, is defined by the expression (see [[#References|[10]]])
+
\frac{h ( x _ {k} ) }{f ( x _ {k} ) }
 +
.
 +
$$
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487069.png" /></td> </tr></table>
+
Simultaneously it is possible to estimate the mean-square error in  $  J _ {N} $,
 +
that is, the quantity  $  \sigma _ {N} = ( {\mathsf D} \zeta / N )  ^ {1/2} $,
 +
and to approximately construct a suitable confidence interval for  $  J $.  
 +
By the choice of the density  $  f $
 +
it is possible to arrange for estimates with possibly smaller variance. For example, if  $  0 \leq  m _ {1} \leq  h / f \leq  m _ {2} < + \infty $,
 +
then  $  {\mathsf D} \zeta \leq  ( m _ {2} - m _ {1} )  ^ {2} / 4 $
 +
and if  $  f = h / J $,
 +
then  $  {\mathsf D} \zeta = 0 $.
 +
The corresponding algorithm is called essential sampling (choice by importance). Another common modification — the method of selection of principal part — occurs when a function  $  h _ {0} \approx h $
 +
is determined whose integral is known. It is sometimes useful to combine Monte-Carlo methods and classical quadratures (cf. [[Quadrature|Quadrature]]) in so-called random quadrature formulas, the basic idea of which is that the nodes and coefficients in any quadrature sum (for example, in interpolation) are chosen randomly from a distribution, providing an [[Unbiased estimator|unbiased estimator]] of the integral [[#References|[3]]]. Particular cases of these formulas are: the so-called method of stratified sampling, in which the nodes are chosen one in each part of a fixed partition of the domain of integration and the coefficients are proportional to the corresponding volumes; and the so-called method of symmetric sampling, which, in the case of integration over  $  ( 0 , 1 ) $,
 +
is defined by the expression (see [[#References|[10]]])
 +
 
 +
$$
 +
2 J  =  {\mathsf E}
 +
\left [
 +
 
 +
\frac{h ( \xi ) + h ( 1 - \xi ) }{f _  \xi  ( \xi ) }
 +
 
 +
\right ] .
 +
$$
  
 
Here the order of the speed of convergence of the Monte-Carlo method is increased and, in certain cases, becomes best possible in the class of problems being considered.
 
Here the order of the speed of convergence of the Monte-Carlo method is increased and, in certain cases, becomes best possible in the class of problems being considered.
Line 59: Line 174:
 
A number of modifications of Monte-Carlo methods are based on the (perhaps formal) representation of the required value as a double integral
 
A number of modifications of Monte-Carlo methods are based on the (perhaps formal) representation of the required value as a double integral
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487070.png" /></td> </tr></table>
+
$$
 +
J  =  \int\limits _ { X } \int\limits _ { Y }
 +
f ( x , y ) h ( x , y )  d x  d y  = \
 +
{\mathsf E} \zeta ,
 +
$$
 +
 
 +
where  $  \zeta = h ( \xi , \eta ) $
 +
and the vector  $  ( \xi , \eta ) $
 +
is distributed with density  $  f ( x , y ) $.
 +
It is known that  $  {\mathsf E} ( \zeta ) = {\mathsf E} {\mathsf E} ( \zeta \mid  \xi ) $
 +
and that
 +
 
 +
$$ \tag{1 }
 +
{\mathsf D} \zeta  = {\mathsf D} {\mathsf E}
 +
( \zeta  \mid  \xi ) + {\mathsf E} {\mathsf D}
 +
( \zeta  \mid  \xi )  = A _ {1} + A _ {2} ,
 +
$$
 +
 
 +
where  $  {\mathsf E} ( \zeta \mid  \xi ) $
 +
is the conditional mathematical expectation and  $  {\mathsf D} ( \zeta \mid  \xi ) $
 +
is the conditional variance of  $  \zeta $
 +
given a fixed value of  $  \xi $.
 +
Formula (1) is widely used in Monte-Carlo methods. In particular, it shows that  $  {\mathsf D} {\mathsf E} ( \zeta \mid  \xi ) < {\mathsf D} \zeta $,
 +
that is, analytic averaging over any variable increases the accuracy of the Monte-Carlo method. However, in this connection, the amount of computation may be significantly increased. The computer time necessary to achieve a given accuracy is proportional to  $  t {\mathsf D} \zeta $,
 +
where  $  t $
 +
is the average time it takes to obtain one value of  $  \zeta $.  
 +
The method of splitting is optimal with respect to this criterion. Its simplest version uses the  "unbiased" estimator
 +
 
 +
$$
 +
\zeta _ {n}  = \
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487071.png" /> and the vector <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487072.png" /> is distributed with density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487073.png" />. It is known that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487074.png" /> and that
+
\frac{1}{n}
 +
\sum _ { k= } 1 ^ { n }  h ( \xi , \eta _ {k} ) ,
 +
$$
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487075.png" /></td> <td valign="top" style="width:5%;text-align:right;">(1)</td></tr></table>
+
where  $  \eta _ {1} \dots \eta _ {n} $
 +
are conditionally-independent and distributed as  $  \eta $
 +
for a fixed value of  $  \xi $.  
 +
Using (1) it is possible to obtain an optimal value
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487076.png" /> is the conditional mathematical expectation and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487077.png" /> is the conditional variance of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487078.png" /> given a fixed value of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487079.png" />. Formula (1) is widely used in Monte-Carlo methods. In particular, it shows that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487080.png" />, that is, analytic averaging over any variable increases the accuracy of the Monte-Carlo method. However, in this connection, the amount of computation may be significantly increased. The computer time necessary to achieve a given accuracy is proportional to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487081.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487082.png" /> is the average time it takes to obtain one value of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487083.png" />. The method of splitting is optimal with respect to this criterion. Its simplest version uses the  "unbiased" estimator
+
$$
 +
\left [
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487084.png" /></td> </tr></table>
+
\frac{A _ {2} }{A _ {1} }
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487085.png" /> are conditionally-independent and distributed as <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487086.png" /> for a fixed value of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487087.png" />. Using (1) it is possible to obtain an optimal value
+
\frac{t _ {1} }{t _ {2} }
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487088.png" /></td> </tr></table>
+
\right ]  ^ {1/2} ,
 +
$$
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487089.png" /> are the average computing times corresponding to the samples <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487090.png" /> (see, for example, [[#References|[4]]]).
+
where $  t _ {1} , t _ {2} $
 +
are the average computing times corresponding to the samples $  \xi , \eta $(
 +
see, for example, [[#References|[4]]]).
  
If the integrand depends on a parameter, it is expedient to use the method of dependent trials, that is, to estimate the integral for various values of the parameter at one and the same random points [[#References|[20]]]. An important property of the Monte-Carlo method is the comparatively relatively weak dependence of the mean-square error <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487091.png" /> on the number of measurements; moreover, the order of convergence relative to the number of points <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487092.png" /> is always one and the same: <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487093.png" />. This allows one to estimate (after a preliminary transformation of the problem) integrals of very high, and even infinite, multiplicity. For example, methods have been worked out for the estimation of Wiener integrals [[#References|[19]]].
+
If the integrand depends on a parameter, it is expedient to use the method of dependent trials, that is, to estimate the integral for various values of the parameter at one and the same random points [[#References|[20]]]. An important property of the Monte-Carlo method is the comparatively relatively weak dependence of the mean-square error $  \sigma _ {N} $
 +
on the number of measurements; moreover, the order of convergence relative to the number of points $  N $
 +
is always one and the same: $  N  ^ {-} 1/2 $.  
 +
This allows one to estimate (after a preliminary transformation of the problem) integrals of very high, and even infinite, multiplicity. For example, methods have been worked out for the estimation of Wiener integrals [[#References|[19]]].
  
 
==Monte-Carlo algorithms for solving integral equations of the second kind.==
 
==Monte-Carlo algorithms for solving integral equations of the second kind.==
Suppose it is required to estimate a linear functional <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487094.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487095.png" />, where the integral operator <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487096.png" /> with kernel <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487097.png" /> satisfies a condition providing convergence of the [[Neumann series|Neumann series]]: <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487098.png" />. A [[Markov chain|Markov chain]] <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m06487099.png" /> is defined by an initial density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870100.png" /> and a transition density <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870101.png" />; the probability of termination of the chain at <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870102.png" /> is equal to
+
Suppose it is required to estimate a linear functional $  J _ {h} = ( \phi , h ) $,  
 +
where $  \phi = K \phi + f $,  
 +
where the integral operator $  K $
 +
with kernel $  k ( x  ^  \prime  , x ) $
 +
satisfies a condition providing convergence of the [[Neumann series|Neumann series]]: $  \| K ^ {n _ {0} } \| < 1 $.  
 +
A [[Markov chain|Markov chain]] $  \{ x _ {n} \} $
 +
is defined by an initial density $  \pi ( x) $
 +
and a transition density $  p ( x  ^  \prime  , x ) = p ( x  ^  \prime  \rightarrow x ) $;  
 +
the probability of termination of the chain at $  x  ^  \prime  $
 +
is equal to
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870103.png" /></td> </tr></table>
+
$$
 +
g ( x  ^  \prime  )  = 1 -
 +
\int\limits p ( x  ^  \prime  , x )  d x .
 +
$$
  
Let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870104.png" /> be the random index of the last state. Further, let a functional of the trajectories of the chain with expectation <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870105.png" /> be defined. Most often the so-called collision estimator
+
Let $  N $
 +
be the random index of the last state. Further, let a functional of the trajectories of the chain with expectation $  J _ {h} $
 +
be defined. Most often the so-called collision estimator
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870106.png" /></td> </tr></table>
+
$$
 +
\xi  = \sum _ { n= } 0 ^ { N }
 +
Q _ {n} h ( x _ {n} )
 +
$$
  
 
is used, where
 
is used, where
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870107.png" /></td> </tr></table>
+
$$
 +
Q _ {0}  = \
 +
 
 +
\frac{f ( x _ {0} ) }{\pi ( x _ {0} ) }
 +
,\ \
 +
Q _ {n}  = Q _ {n-} 1
  
If <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870108.png" /> for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870109.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870110.png" /> for <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870111.png" />, then under certain additional conditions
+
\frac{k ( x _ {n-} 1 , x _ {n} ) }{p ( x _ {n-} 1 , x _ {n} ) }
 +
.
 +
$$
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870112.png" /></td> </tr></table>
+
If  $  p ( x  ^  \prime  , x ) \neq 0 $
 +
for  $  k ( x  ^  \prime  , x ) \neq 0 $
 +
and  $  \pi ( x) \neq 0 $
 +
for  $  f ( x) \neq 0 $,
 +
then under certain additional conditions
 +
 
 +
$$
 +
{\mathsf E} \xi  = \
 +
\sum _ { n= } 0 ^  \infty 
 +
( K  ^ {n} f , h )  = \
 +
( \phi , h )  = \
 +
\int\limits _ { X } \phi ( x) h ( x)  d x
 +
$$
  
 
(see [[#References|[3]]]–[[#References|[5]]]). The possibility of attaining a small variance in the case of constant sign is ensured by the following result: If
 
(see [[#References|[3]]]–[[#References|[5]]]). The possibility of attaining a small variance in the case of constant sign is ensured by the following result: If
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870113.png" /></td> </tr></table>
+
$$
 +
\pi ( x)  = \
 +
 
 +
\frac{f ( x) \phi  ^ {*} ( x) }{( f , \phi  ^ {*} ) }
 +
\ \
 +
\textrm{ and } \ \
 +
p ( x  ^  \prime  , x )  = \
 +
 
 +
\frac{k ( x  ^  \prime  , x ) \phi  ^ {*} ( x) }{[ K  ^ {*} \phi  ^ {*} ] ( x  ^  \prime  ) }
 +
,
 +
$$
  
where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870114.png" />, then <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870115.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870116.png" /> (see [[#References|[4]]]). If a suitable Markov chain is simulated on a computer, statistical estimates of linear functionals in the solution of an integral equation of the second kind can be obtained. This gives the possibility of locally estimating the solution on the basis of the representation <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870117.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870118.png" />. In a number of cases, alongside Monte-Carlo methods, number-theoretic methods are applied in order to solve these problems (see [[#References|[21]]]). A Monte-Carlo method for estimating the first eigen value of an integral operator has been realized by an iteration method on the basis of the relation [[#References|[22]]]
+
where $  \phi  ^ {*} = K  ^ {*} \phi  ^ {*} + h $,  
 +
then $  {\mathsf D} \xi = 0 $
 +
and $  {\mathsf E} \xi = J _ {h} $(
 +
see [[#References|[4]]]). If a suitable Markov chain is simulated on a computer, statistical estimates of linear functionals in the solution of an integral equation of the second kind can be obtained. This gives the possibility of locally estimating the solution on the basis of the representation $  \phi ( x) = ( \phi , h _ {x} ) + f ( x) $,  
 +
where $  h _ {x} ( x  ^  \prime  ) = k ( x  ^  \prime  , x ) $.  
 +
In a number of cases, alongside Monte-Carlo methods, number-theoretic methods are applied in order to solve these problems (see [[#References|[21]]]). A Monte-Carlo method for estimating the first eigen value of an integral operator has been realized by an iteration method on the basis of the relation [[#References|[22]]]
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870119.png" /></td> </tr></table>
+
$$
 +
{\mathsf E} [ Q _ {n} h ( x _ {n} ) ]  = \
 +
( K  ^ {n} f , h ) .
 +
$$
  
All these results can be almost automatically extended to systems of linear algebraic equations of the form <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870120.png" /> (see [[#References|[23]]]).
+
All these results can be almost automatically extended to systems of linear algebraic equations of the form $  x + H x = h $(
 +
see [[#References|[23]]]).
  
 
==Modifications of Monte-Carlo methods in radiative transport theory.==
 
==Modifications of Monte-Carlo methods in radiative transport theory.==
(See [[#References|[11]]]–[[#References|[17]]].) The density of the average number of particle collisions in a phase space with coordinates <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870121.png" /> and velocities <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870122.png" /> satisfies an integral equation of the second kind; its kernel in the single-velocity case has the form
+
(See [[#References|[11]]]–[[#References|[17]]].) The density of the average number of particle collisions in a phase space with coordinates $  \mathbf r $
 +
and velocities $  \omega $
 +
satisfies an integral equation of the second kind; its kernel in the single-velocity case has the form
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870123.png" /></td> </tr></table>
+
$$
  
Here <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870124.png" /> is the scattering coefficient (section), <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870125.png" /> is the relaxation coefficient, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870126.png" /> is the scattering indicatrix, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870127.png" /> is the optical length of a path from <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870128.png" /> to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870129.png" /> (see [[#References|[3]]], [[#References|[4]]]). For the construction of estimates with small variances one uses, for example, asymptotic solutions of the adjoint transport equation [[#References|[4]]]; the simplest algorithm of this type is the so-called exponential transformation (see [[#References|[4]]], [[#References|[11]]]). A modification of a local estimate of the flow of particles has been developed (see [[#References|[3]]], [[#References|[4]]], [[#References|[11]]]–[[#References|[13]]], [[#References|[17]]], [[#References|[18]]]). Using a simulation of a Markov chain (for example, the physical transport process in a certain medium) it is possible to simultaneously obtain dependent estimators of functionals for different values of the parameters; by differentiating the  "weight"  <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870130.png" /> it is sometimes possible to obtain unbiased estimators of the corresponding derivatives (see [[#References|[4]]], [[#References|[12]]]). This provides an opportunity to use Monte-Carlo methods in solving certain inverse problems [[#References|[12]]]. For the solution of problems in transport theory  "bifurcation"  of the trajectory and analytic averaging is effective [[#References|[11]]]. The simulation of trajectories of particles in compound media can sometimes be essentially simplified by the method of a maximal section (see [[#References|[3]]]–[[#References|[5]]]).
+
\frac{\sigma _ {s} ( \mathbf r  ^  \prime  ) g ( \mu )
 +
\mathop{\rm exp} ( - \tau ( \mathbf r  ^  \prime  , \mathbf r ) ) \sigma (
 +
\mathbf r ) }{\sigma ( \mathbf r  ^  \prime  ) 2 \pi
 +
| \mathbf r  ^  \prime  - \mathbf r |  ^ {2} }
 +
 
 +
\delta
 +
\left (
 +
\omega -
 +
 
 +
\frac{\mathbf r - \mathbf r  ^  \prime  }{| \mathbf r - \mathbf r  ^  \prime  | }
 +
 
 +
\right ) .
 +
$$
 +
 
 +
Here $  \sigma _ {s} ( \mathbf r ) $
 +
is the scattering coefficient (section), $  \sigma ( \mathbf r ) $
 +
is the relaxation coefficient, $  g ( \mu ) $
 +
is the scattering indicatrix, and $  \tau ( \mathbf r  ^  \prime  , \mathbf r ) $
 +
is the optical length of a path from $  \mathbf r  ^  \prime  $
 +
to $  \mathbf r $(
 +
see [[#References|[3]]], [[#References|[4]]]). For the construction of estimates with small variances one uses, for example, asymptotic solutions of the adjoint transport equation [[#References|[4]]]; the simplest algorithm of this type is the so-called exponential transformation (see [[#References|[4]]], [[#References|[11]]]). A modification of a local estimate of the flow of particles has been developed (see [[#References|[3]]], [[#References|[4]]], [[#References|[11]]]–[[#References|[13]]], [[#References|[17]]], [[#References|[18]]]). Using a simulation of a Markov chain (for example, the physical transport process in a certain medium) it is possible to simultaneously obtain dependent estimators of functionals for different values of the parameters; by differentiating the  "weight"   $ Q _ {n} $
 +
it is sometimes possible to obtain unbiased estimators of the corresponding derivatives (see [[#References|[4]]], [[#References|[12]]]). This provides an opportunity to use Monte-Carlo methods in solving certain inverse problems [[#References|[12]]]. For the solution of problems in transport theory  "bifurcation"  of the trajectory and analytic averaging is effective [[#References|[11]]]. The simulation of trajectories of particles in compound media can sometimes be essentially simplified by the method of a maximal section (see [[#References|[3]]]–[[#References|[5]]]).
  
 
These are constructed on the basis of the corresponding integral relations. For example, the standard five-point difference approximation for the Laplace equation has the form of the formula for the complete mathematical expectation of a symmetric random walk over a grid with absorption at the boundary (see, for example, [[#References|[2]]], [[#References|[3]]]). A continuous analogue of this formula is the relation
 
These are constructed on the basis of the corresponding integral relations. For example, the standard five-point difference approximation for the Laplace equation has the form of the formula for the complete mathematical expectation of a symmetric random walk over a grid with absorption at the boundary (see, for example, [[#References|[2]]], [[#References|[3]]]). A continuous analogue of this formula is the relation
  
<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870131.png" /></td> <td valign="top" style="width:5%;text-align:right;">(2)</td></tr></table>
+
$$ \tag{2 }
 +
u ( P) = \
  
where the integral is taken over the surface of a sphere lying entirely within the given domain and with centre at <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/m/m064/m064870/m064870132.png" />. Formula (2) and other similar relations provide an opportunity to use isotropic  "random walk on spheres"  when solving elliptic and parabolic equations (see [[#References|[24]]], [[#References|[4]]]). Monte-Carlo methods are effective, for example, for estimating the solution of multi-dimensional boundary value problems at a point.
+
\frac{\int\limits _ {S ( P) } u ( \mathbf r ( s) )  d s }{\int\limits _ {S ( P) } d s }
 +
,
 +
$$
 +
 
 +
where the integral is taken over the surface of a sphere lying entirely within the given domain and with centre at $  P $.  
 +
Formula (2) and other similar relations provide an opportunity to use isotropic  "random walk on spheres"  when solving elliptic and parabolic equations (see [[#References|[24]]], [[#References|[4]]]). Monte-Carlo methods are effective, for example, for estimating the solution of multi-dimensional boundary value problems at a point.
  
 
A simulation of Markov branching processes allows one to construct estimates of the solution of certain non-linear equations, for example, the Boltzmann equation in the theory of rarefied gases [[#References|[3]]].
 
A simulation of Markov branching processes allows one to construct estimates of the solution of certain non-linear equations, for example, the Boltzmann equation in the theory of rarefied gases [[#References|[3]]].
Line 121: Line 361:
 
====References====
 
====References====
 
<table><TR><TD valign="top">[1]</TD> <TD valign="top">  J. von Neumann,  ''Nat. Bureau Standard Appl. Math. Series'' , '''12'''  (1951)  pp. 36–38</TD></TR><TR><TD valign="top">[2]</TD> <TD valign="top">  N.P. Buslenko,  et al.,  "The method of statistical trials (the Monte-Carlo method)" , Moscow  (1962)  (In Russian)</TD></TR><TR><TD valign="top">[3]</TD> <TD valign="top">  S.M. Ermakov,  "Die Monte-Carlo Methode und verwandte Fragen" , Deutsch. Verlag Wissenschaft.  (1975)  (Translated from Russian)</TD></TR><TR><TD valign="top">[4]</TD> <TD valign="top">  G.A. Mikhailov,  "Some questions in the theory of Monte-Carlo methods" , Novosibirsk  (1971)  (In Russian)</TD></TR><TR><TD valign="top">[5]</TD> <TD valign="top">  I.M. Sobol',  "Numerical Monte-Carlo methods" , Moscow  (1973)  (In Russian)</TD></TR><TR><TD valign="top">[6]</TD> <TD valign="top">  Yu.G. Pollyak,  "Probabilistic simulation on computers" , Moscow  (1971)  (In Russian)</TD></TR><TR><TD valign="top">[7]</TD> <TD valign="top">  N.S. Bakhvalov,  "On optimal convergence estimates for quadrature processes and integration methods of Monte-Carlo type on function classes" , ''Numerical Methods for Solving Differential and Integral Equations and Quadrature Formulas'' , Moscow  (1964)  pp. 5–63  (In Russian)</TD></TR><TR><TD valign="top">[8]</TD> <TD valign="top">  N.S. Bakhvalov,  "An estimate of the remainder term in quadrature formula"  ''Comp. Math. Math. Phys.'' , '''1''' :  1  (1961)  pp. 68–82  ''Zh. Vychisl. Mat. i Mat. Fiz.'' , '''1''' :  1  (1961)  pp. 64–77</TD></TR><TR><TD valign="top">[9]</TD> <TD valign="top">  N.S. Bakhvalov,  "Approximate computation of multiple integrals"  ''Vestnik Moskov. Gos. Univ. Ser. Mat. Mekh. Astronom. Fiz. Khim.'' , '''4'''  (1959)  pp. 3–18  (In Russian)</TD></TR><TR><TD valign="top">[10]</TD> <TD valign="top">  J.M. Hammersley,  D.C. Handscomb,  "Monte-Carlo methods" , Methuen  (1964)</TD></TR><TR><TD valign="top">[11]</TD> <TD valign="top"> , ''Monte-Carlo methods and problems of radiative transport'' , Moscow  (1967)</TD></TR><TR><TD valign="top">[12]</TD> <TD valign="top">  G.I. Marchuk,  et al.,  "The Monte-Carlo method in atmospheric optics" , Novosibirsk  (1976)  (In Russian)</TD></TR><TR><TD valign="top">[13]</TD> <TD valign="top">  J. Spanier,  E. Gelbard,  "Monte-Carlo principles and neutron transport problems" , Addison-Wesley  (1969)</TD></TR><TR><TD valign="top">[14]</TD> <TD valign="top">  V.V. Chavchanidze,  ''Izv. Akad. Nauk SSSR Ser. Fiz.'' , '''19''' :  6  (1955)  pp. 629–638</TD></TR><TR><TD valign="top">[15]</TD> <TD valign="top"> , ''The penetration of radiation through non-uniformity in protection'' , Moscow  (1968)  (In Russian)</TD></TR><TR><TD valign="top">[16]</TD> <TD valign="top">  A.D. Frank-Kamenetskii,  ''Atomnaya Energiya'' , '''16''' :  2  (1964)  pp. 119–122</TD></TR><TR><TD valign="top">[17]</TD> <TD valign="top">  M.H. Kalos,  ''Nuclear Sci. and Eng.'' , '''33'''  (1968)  pp. 284–290</TD></TR><TR><TD valign="top">[18]</TD> <TD valign="top"> , ''Monte-Carlo methods and their application. Reports III All-Union Conf. Monte-Carlo Methods'' , Novosibirsk  (1971)</TD></TR><TR><TD valign="top">[19]</TD> <TD valign="top">  I.M. Gelfand,  A.S. Frolov,  N.N. Chentsov,  "The computation of continuous integrals by the Monte-Carlo method"  ''Izv. Vuz. Ser. Mat.'' , '''5'''  (1958)  pp. 32–45  (In Russian)</TD></TR><TR><TD valign="top">[20]</TD> <TD valign="top">  A.S. Frolov,  N.N. Chentsov,  "On the calculation of definite integrals dependent on a parameter by the Monte-Carlo method"  ''USSR Comp. Math. Math. Phys.'' , '''2''' :  4  (1962)  pp. 802–807  ''Zh. Vychisl. Mat. i. Mat. Fiz.'' , '''2''' :  4  (1962)  pp. 714–717</TD></TR><TR><TD valign="top">[21]</TD> <TD valign="top">  N.M. Korobov,  "Number-theoretic methods in applied analysis" , Moscow  (1963)  (In Russian)</TD></TR><TR><TD valign="top">[22]</TD> <TD valign="top">  V.S. Vladimirov,  "Monte-Carlo methods as applied to the calculation of the lowest eigenvalue and the associated eigenfunction of a linear integral equation"  ''Theor. Probab. Appl.'' , '''1''' :  1  (1956)  pp. 101–116  ''Teor. Veroyatnost. i. Primenen.'' , '''1''' :  1  (1956)  pp. 113–130</TD></TR><TR><TD valign="top">[23]</TD> <TD valign="top">  J.H. Curtiss,  "Monte-Carlo methods for the iteration of linear operators"  ''J. Math. Phys.'' , '''32''' :  4  (1954)  pp. 209–232</TD></TR><TR><TD valign="top">[24]</TD> <TD valign="top">  M.E. Muller,  "Some continuous Monte-Carlo methods for the Dirichlet problem"  ''Ann. Math. Stat.'' , '''27''' :  3  (1956)  pp. 569–589</TD></TR></table>
 
<table><TR><TD valign="top">[1]</TD> <TD valign="top">  J. von Neumann,  ''Nat. Bureau Standard Appl. Math. Series'' , '''12'''  (1951)  pp. 36–38</TD></TR><TR><TD valign="top">[2]</TD> <TD valign="top">  N.P. Buslenko,  et al.,  "The method of statistical trials (the Monte-Carlo method)" , Moscow  (1962)  (In Russian)</TD></TR><TR><TD valign="top">[3]</TD> <TD valign="top">  S.M. Ermakov,  "Die Monte-Carlo Methode und verwandte Fragen" , Deutsch. Verlag Wissenschaft.  (1975)  (Translated from Russian)</TD></TR><TR><TD valign="top">[4]</TD> <TD valign="top">  G.A. Mikhailov,  "Some questions in the theory of Monte-Carlo methods" , Novosibirsk  (1971)  (In Russian)</TD></TR><TR><TD valign="top">[5]</TD> <TD valign="top">  I.M. Sobol',  "Numerical Monte-Carlo methods" , Moscow  (1973)  (In Russian)</TD></TR><TR><TD valign="top">[6]</TD> <TD valign="top">  Yu.G. Pollyak,  "Probabilistic simulation on computers" , Moscow  (1971)  (In Russian)</TD></TR><TR><TD valign="top">[7]</TD> <TD valign="top">  N.S. Bakhvalov,  "On optimal convergence estimates for quadrature processes and integration methods of Monte-Carlo type on function classes" , ''Numerical Methods for Solving Differential and Integral Equations and Quadrature Formulas'' , Moscow  (1964)  pp. 5–63  (In Russian)</TD></TR><TR><TD valign="top">[8]</TD> <TD valign="top">  N.S. Bakhvalov,  "An estimate of the remainder term in quadrature formula"  ''Comp. Math. Math. Phys.'' , '''1''' :  1  (1961)  pp. 68–82  ''Zh. Vychisl. Mat. i Mat. Fiz.'' , '''1''' :  1  (1961)  pp. 64–77</TD></TR><TR><TD valign="top">[9]</TD> <TD valign="top">  N.S. Bakhvalov,  "Approximate computation of multiple integrals"  ''Vestnik Moskov. Gos. Univ. Ser. Mat. Mekh. Astronom. Fiz. Khim.'' , '''4'''  (1959)  pp. 3–18  (In Russian)</TD></TR><TR><TD valign="top">[10]</TD> <TD valign="top">  J.M. Hammersley,  D.C. Handscomb,  "Monte-Carlo methods" , Methuen  (1964)</TD></TR><TR><TD valign="top">[11]</TD> <TD valign="top"> , ''Monte-Carlo methods and problems of radiative transport'' , Moscow  (1967)</TD></TR><TR><TD valign="top">[12]</TD> <TD valign="top">  G.I. Marchuk,  et al.,  "The Monte-Carlo method in atmospheric optics" , Novosibirsk  (1976)  (In Russian)</TD></TR><TR><TD valign="top">[13]</TD> <TD valign="top">  J. Spanier,  E. Gelbard,  "Monte-Carlo principles and neutron transport problems" , Addison-Wesley  (1969)</TD></TR><TR><TD valign="top">[14]</TD> <TD valign="top">  V.V. Chavchanidze,  ''Izv. Akad. Nauk SSSR Ser. Fiz.'' , '''19''' :  6  (1955)  pp. 629–638</TD></TR><TR><TD valign="top">[15]</TD> <TD valign="top"> , ''The penetration of radiation through non-uniformity in protection'' , Moscow  (1968)  (In Russian)</TD></TR><TR><TD valign="top">[16]</TD> <TD valign="top">  A.D. Frank-Kamenetskii,  ''Atomnaya Energiya'' , '''16''' :  2  (1964)  pp. 119–122</TD></TR><TR><TD valign="top">[17]</TD> <TD valign="top">  M.H. Kalos,  ''Nuclear Sci. and Eng.'' , '''33'''  (1968)  pp. 284–290</TD></TR><TR><TD valign="top">[18]</TD> <TD valign="top"> , ''Monte-Carlo methods and their application. Reports III All-Union Conf. Monte-Carlo Methods'' , Novosibirsk  (1971)</TD></TR><TR><TD valign="top">[19]</TD> <TD valign="top">  I.M. Gelfand,  A.S. Frolov,  N.N. Chentsov,  "The computation of continuous integrals by the Monte-Carlo method"  ''Izv. Vuz. Ser. Mat.'' , '''5'''  (1958)  pp. 32–45  (In Russian)</TD></TR><TR><TD valign="top">[20]</TD> <TD valign="top">  A.S. Frolov,  N.N. Chentsov,  "On the calculation of definite integrals dependent on a parameter by the Monte-Carlo method"  ''USSR Comp. Math. Math. Phys.'' , '''2''' :  4  (1962)  pp. 802–807  ''Zh. Vychisl. Mat. i. Mat. Fiz.'' , '''2''' :  4  (1962)  pp. 714–717</TD></TR><TR><TD valign="top">[21]</TD> <TD valign="top">  N.M. Korobov,  "Number-theoretic methods in applied analysis" , Moscow  (1963)  (In Russian)</TD></TR><TR><TD valign="top">[22]</TD> <TD valign="top">  V.S. Vladimirov,  "Monte-Carlo methods as applied to the calculation of the lowest eigenvalue and the associated eigenfunction of a linear integral equation"  ''Theor. Probab. Appl.'' , '''1''' :  1  (1956)  pp. 101–116  ''Teor. Veroyatnost. i. Primenen.'' , '''1''' :  1  (1956)  pp. 113–130</TD></TR><TR><TD valign="top">[23]</TD> <TD valign="top">  J.H. Curtiss,  "Monte-Carlo methods for the iteration of linear operators"  ''J. Math. Phys.'' , '''32''' :  4  (1954)  pp. 209–232</TD></TR><TR><TD valign="top">[24]</TD> <TD valign="top">  M.E. Muller,  "Some continuous Monte-Carlo methods for the Dirichlet problem"  ''Ann. Math. Stat.'' , '''27''' :  3  (1956)  pp. 569–589</TD></TR></table>
 
 
  
 
====Comments====
 
====Comments====

Latest revision as of 08:01, 6 June 2020


method of statistical trials

A numerical method based on simulation by random variables and the construction of statistical estimators for the unknown quantities. It is usually supposed that the Monte-Carlo method originated in 1949 (see [1]) when, in connection with work on the construction of atomic reactors, J. von Neumann and S. Ulam suggested using the apparatus of probability theory in the computer solution of applied problems. The Monte-Carlo method is named after the town of Monte-Carlo, famous for its casinos.

Simulation by random variables with given distributions.

As a rule such a simulation is achieved by a transformation of one or more independent values of a random number $ \alpha $, uniformly distributed in the interval $ ( 0 , 1 ) $. The sequence of "sample" values of $ \alpha $ is usually obtained on a computer by number-theoretic algorithms, of which the most widely used is the so-called method of residues, for example, in the form

$$ u _ {0} = 1 ,\ \ u _ {n} \equiv u _ {n-} 1 5 ^ {2p+} 1 \ ( \mathop{\rm mod} 2 ^ {m} ) ,\ \ \alpha _ {n} = u _ {n} \cdot 2 ^ {-} m . $$

Here $ m $ is the order of the mantissa of the computer and

$$ p = \max \{ {q } : {5 ^ {2q+} 1 < 2 ^ {m} } \} . $$

Numbers of this type are called pseudo-random numbers; they are used in statistical testing and in solving typical problems (see [2][6]). The length of the period of the above version of the method of residues is $ 2 ^ {m-} 2 $. Physical generators, tables of random numbers and quasi-random numbers are also used in the Monte-Carlo method. There are Monte-Carlo methods with a small number of playing parameters (see [7]).

The standard method for simulating a discrete random variable with distribution $ {\mathsf P} \{ \xi = x _ {k} \} = p _ {k} $, $ k = 0 , 1 \dots $ is as follows: Put $ \xi = x _ {m} $ if the chosen value of $ \alpha $ satisfies

$$ \sum _ { k= } 0 ^ { m- } 1 p _ {k} \leq \alpha < \ \sum _ { k= } 0 ^ { m } p _ {k} . $$

The standard method for simulating a continuous random variable (sometimes called the method of inverse functions) is to use the, easily verified, representation $ \xi = F ^ { - 1 } ( \alpha ) $, where $ F $ is the distribution function with given density $ f $. Sometimes randomization of the simulation is useful (in other words, the method of superposition), based on the expression

$$ f ( x) = \sum _ { k } p _ {k} f _ {k} ( x) ; $$

here one first chooses a number $ m $ with distribution $ {\mathsf P} \{ m = k \} = p _ {k} $, and then obtains a sample value $ \xi $ from the distribution with density $ f _ {m} $. In other methods of randomization certain parameters of a deterministic method for solving the problem are considered as random variables (see [7][9]).

Another, more general, method for simulating a continuous random variable is the method of exclusion (method of selection), at the basis of which lies the following result: If $ ( \xi , \eta ) $ is uniformly distributed in a domain $ G = \{ {( x , y ) } : {0 \leq y \leq g ( x) } \} $, then $ f _ \xi ( x) = g ( x) / \overline{G}\; $. In the method of exclusion, choose a point $ ( \xi _ {0} , \eta ) $ uniformly in a domain $ G _ {1} \supset G $ and put $ \xi = \xi _ {0} $ if $ ( \xi _ {0} , \eta ) \in G $; otherwise repeat the selection of $ ( \xi _ {0} , \eta ) $, etc. For example, if $ a \leq \xi \leq b $ and $ g ( x) = c f ( x) \leq R $, one can take $ \xi _ {0} = a + ( b - a ) \alpha _ {1} $, $ \eta = R \alpha _ {2} $. The average number of operations in the method of exclusion is proportional to $ \overline{G}\; _ {1} / \overline{G}\; $.

For many random variables, special representations of the form $ \xi = \phi ( \alpha _ {1} \dots \alpha _ {n} ) $ have been obtained. For example, the random variables

$$ \sqrt {- 2 \mathop{\rm ln} \alpha _ {1} } \cdot \cos 2 \pi \alpha _ {2} \ \ \textrm{ and } \ \ \sqrt {- 2 \mathop{\rm ln} \alpha _ {2} } \cdot \sin 2 \pi \alpha _ {2} $$

have standard normal distributions and are independent; the random variable $ \mathop{\rm ln} ( \alpha _ {1} \dots \alpha _ {n} ) $ has the gamma-distribution with parameter $ n $; the random variable $ \max ( \alpha _ {1} \dots \alpha _ {n} ) $ is distributed with density $ n x ^ {n-} 1 $, $ 0 \leq x \leq 1 $; the random variable $ \mathop{\rm exp} \{ \sum _ {k=} 1 ^ {n} ( \mathop{\rm ln} \alpha ) / ( p + k + 1 ) \} $ has the beta-distribution with parameters $ p , n $( see [3][6]).

The standard algorithm for simulating a continuous random vector $ \\vec{xi} = ( \xi _ {1} \dots \xi _ {n} ) $ is to successively choose the values of its components from conditional distributions corresponding to the representation

$$ f {} _ {\\vec{xi} } ( x _ {1} \dots x _ {n} ) = \ f _ {1} ( x _ {1} ) f _ {2} ( x _ {2} \mid x _ {1} ) \dots f _ {n} ( x _ {n} \mid x _ {1} \dots x _ {n-} 1 ) . $$

The method of exclusion extends to the multi-dimensional case without change; it is only necessary, in its formulation, to regard $ \xi $, $ \xi _ {0} $ and $ x $ as vectors. A multi-dimensional normal vector can be simulated by using a special linear transformation of a vector of independent standard normal random variables. Special methods have also been developed for the approximate simulation of stationary Gaussian processes (see, for example, [3], [6]).

If, in a Monte-Carlo method calculation, random variables defined by a real phenomenon are simulated, then the calculation is a direct simulation (imitation) of this phenomenon. Computer simulations have been worked out for: the processes of transport, scattering and reproduction of particles: neutrons, gamma-quanta, photons, electrons, and others (see, for example, [11][18]); simulations of the evolution of ensembles of molecules for the solution of various problems in classical and quantum statistical physics (see, for example, [10][18]); simulation of queueing and industrial processes (see, for example [2], [6], [18]); simulation of various random processes in technology, biology; etc. (see [18]). Simulation algorithms are usually carefully developed; for example, they tabulate complicated functions, modify standard procedures, etc. None-the-less, a direct simulation often cannot provide the necessary accuracy in the estimates of the required variables. Many methods have been developed for increasing the effectiveness of a simulation.

Monte-Carlo algorithms for estimating multiple integrals.

Suppose it is required to estimate an integral $ J = \int h ( x) d x $ with respect to the Lebesgue measure in an $ s $- dimensional Euclidean space $ X $ and let $ f _ \xi ( x) $ be a probability density such that $ J $ can be written as:

$$ J = \int\limits _ { X } f _ \xi ( x) \frac{h ( x) }{f _ \xi ( x) } \ d x = {\mathsf E} \zeta , $$

where $ \zeta = h ( \xi ) / f _ \xi ( \xi ) $. By computer simulation of $ \xi $ it is possible to obtain $ N $ sample values $ x _ {1} \dots x _ {N} $. By the law of large numbers,

$$ J \approx J _ {N} = \ \frac{1}{N} \sum _ { k= } 1 ^ { N } \frac{h ( x _ {k} ) }{f ( x _ {k} ) } . $$

Simultaneously it is possible to estimate the mean-square error in $ J _ {N} $, that is, the quantity $ \sigma _ {N} = ( {\mathsf D} \zeta / N ) ^ {1/2} $, and to approximately construct a suitable confidence interval for $ J $. By the choice of the density $ f $ it is possible to arrange for estimates with possibly smaller variance. For example, if $ 0 \leq m _ {1} \leq h / f \leq m _ {2} < + \infty $, then $ {\mathsf D} \zeta \leq ( m _ {2} - m _ {1} ) ^ {2} / 4 $ and if $ f = h / J $, then $ {\mathsf D} \zeta = 0 $. The corresponding algorithm is called essential sampling (choice by importance). Another common modification — the method of selection of principal part — occurs when a function $ h _ {0} \approx h $ is determined whose integral is known. It is sometimes useful to combine Monte-Carlo methods and classical quadratures (cf. Quadrature) in so-called random quadrature formulas, the basic idea of which is that the nodes and coefficients in any quadrature sum (for example, in interpolation) are chosen randomly from a distribution, providing an unbiased estimator of the integral [3]. Particular cases of these formulas are: the so-called method of stratified sampling, in which the nodes are chosen one in each part of a fixed partition of the domain of integration and the coefficients are proportional to the corresponding volumes; and the so-called method of symmetric sampling, which, in the case of integration over $ ( 0 , 1 ) $, is defined by the expression (see [10])

$$ 2 J = {\mathsf E} \left [ \frac{h ( \xi ) + h ( 1 - \xi ) }{f _ \xi ( \xi ) } \right ] . $$

Here the order of the speed of convergence of the Monte-Carlo method is increased and, in certain cases, becomes best possible in the class of problems being considered.

In general, the domain of integration is partitioned into parallelopipeds. In each parallelopiped the value of the integral is calculated as the average of the value at a random point and the point symmetric to it relative to the centre of the parallelopiped.

A number of modifications of Monte-Carlo methods are based on the (perhaps formal) representation of the required value as a double integral

$$ J = \int\limits _ { X } \int\limits _ { Y } f ( x , y ) h ( x , y ) d x d y = \ {\mathsf E} \zeta , $$

where $ \zeta = h ( \xi , \eta ) $ and the vector $ ( \xi , \eta ) $ is distributed with density $ f ( x , y ) $. It is known that $ {\mathsf E} ( \zeta ) = {\mathsf E} {\mathsf E} ( \zeta \mid \xi ) $ and that

$$ \tag{1 } {\mathsf D} \zeta = {\mathsf D} {\mathsf E} ( \zeta \mid \xi ) + {\mathsf E} {\mathsf D} ( \zeta \mid \xi ) = A _ {1} + A _ {2} , $$

where $ {\mathsf E} ( \zeta \mid \xi ) $ is the conditional mathematical expectation and $ {\mathsf D} ( \zeta \mid \xi ) $ is the conditional variance of $ \zeta $ given a fixed value of $ \xi $. Formula (1) is widely used in Monte-Carlo methods. In particular, it shows that $ {\mathsf D} {\mathsf E} ( \zeta \mid \xi ) < {\mathsf D} \zeta $, that is, analytic averaging over any variable increases the accuracy of the Monte-Carlo method. However, in this connection, the amount of computation may be significantly increased. The computer time necessary to achieve a given accuracy is proportional to $ t {\mathsf D} \zeta $, where $ t $ is the average time it takes to obtain one value of $ \zeta $. The method of splitting is optimal with respect to this criterion. Its simplest version uses the "unbiased" estimator

$$ \zeta _ {n} = \ \frac{1}{n} \sum _ { k= } 1 ^ { n } h ( \xi , \eta _ {k} ) , $$

where $ \eta _ {1} \dots \eta _ {n} $ are conditionally-independent and distributed as $ \eta $ for a fixed value of $ \xi $. Using (1) it is possible to obtain an optimal value

$$ n = \left [ \frac{A _ {2} }{A _ {1} } \frac{t _ {1} }{t _ {2} } \right ] ^ {1/2} , $$

where $ t _ {1} , t _ {2} $ are the average computing times corresponding to the samples $ \xi , \eta $( see, for example, [4]).

If the integrand depends on a parameter, it is expedient to use the method of dependent trials, that is, to estimate the integral for various values of the parameter at one and the same random points [20]. An important property of the Monte-Carlo method is the comparatively relatively weak dependence of the mean-square error $ \sigma _ {N} $ on the number of measurements; moreover, the order of convergence relative to the number of points $ N $ is always one and the same: $ N ^ {-} 1/2 $. This allows one to estimate (after a preliminary transformation of the problem) integrals of very high, and even infinite, multiplicity. For example, methods have been worked out for the estimation of Wiener integrals [19].

Monte-Carlo algorithms for solving integral equations of the second kind.

Suppose it is required to estimate a linear functional $ J _ {h} = ( \phi , h ) $, where $ \phi = K \phi + f $, where the integral operator $ K $ with kernel $ k ( x ^ \prime , x ) $ satisfies a condition providing convergence of the Neumann series: $ \| K ^ {n _ {0} } \| < 1 $. A Markov chain $ \{ x _ {n} \} $ is defined by an initial density $ \pi ( x) $ and a transition density $ p ( x ^ \prime , x ) = p ( x ^ \prime \rightarrow x ) $; the probability of termination of the chain at $ x ^ \prime $ is equal to

$$ g ( x ^ \prime ) = 1 - \int\limits p ( x ^ \prime , x ) d x . $$

Let $ N $ be the random index of the last state. Further, let a functional of the trajectories of the chain with expectation $ J _ {h} $ be defined. Most often the so-called collision estimator

$$ \xi = \sum _ { n= } 0 ^ { N } Q _ {n} h ( x _ {n} ) $$

is used, where

$$ Q _ {0} = \ \frac{f ( x _ {0} ) }{\pi ( x _ {0} ) } ,\ \ Q _ {n} = Q _ {n-} 1 \frac{k ( x _ {n-} 1 , x _ {n} ) }{p ( x _ {n-} 1 , x _ {n} ) } . $$

If $ p ( x ^ \prime , x ) \neq 0 $ for $ k ( x ^ \prime , x ) \neq 0 $ and $ \pi ( x) \neq 0 $ for $ f ( x) \neq 0 $, then under certain additional conditions

$$ {\mathsf E} \xi = \ \sum _ { n= } 0 ^ \infty ( K ^ {n} f , h ) = \ ( \phi , h ) = \ \int\limits _ { X } \phi ( x) h ( x) d x $$

(see [3][5]). The possibility of attaining a small variance in the case of constant sign is ensured by the following result: If

$$ \pi ( x) = \ \frac{f ( x) \phi ^ {*} ( x) }{( f , \phi ^ {*} ) } \ \ \textrm{ and } \ \ p ( x ^ \prime , x ) = \ \frac{k ( x ^ \prime , x ) \phi ^ {*} ( x) }{[ K ^ {*} \phi ^ {*} ] ( x ^ \prime ) } , $$

where $ \phi ^ {*} = K ^ {*} \phi ^ {*} + h $, then $ {\mathsf D} \xi = 0 $ and $ {\mathsf E} \xi = J _ {h} $( see [4]). If a suitable Markov chain is simulated on a computer, statistical estimates of linear functionals in the solution of an integral equation of the second kind can be obtained. This gives the possibility of locally estimating the solution on the basis of the representation $ \phi ( x) = ( \phi , h _ {x} ) + f ( x) $, where $ h _ {x} ( x ^ \prime ) = k ( x ^ \prime , x ) $. In a number of cases, alongside Monte-Carlo methods, number-theoretic methods are applied in order to solve these problems (see [21]). A Monte-Carlo method for estimating the first eigen value of an integral operator has been realized by an iteration method on the basis of the relation [22]

$$ {\mathsf E} [ Q _ {n} h ( x _ {n} ) ] = \ ( K ^ {n} f , h ) . $$

All these results can be almost automatically extended to systems of linear algebraic equations of the form $ x + H x = h $( see [23]).

Modifications of Monte-Carlo methods in radiative transport theory.

(See [11][17].) The density of the average number of particle collisions in a phase space with coordinates $ \mathbf r $ and velocities $ \omega $ satisfies an integral equation of the second kind; its kernel in the single-velocity case has the form

$$ \frac{\sigma _ {s} ( \mathbf r ^ \prime ) g ( \mu ) \mathop{\rm exp} ( - \tau ( \mathbf r ^ \prime , \mathbf r ) ) \sigma ( \mathbf r ) }{\sigma ( \mathbf r ^ \prime ) 2 \pi | \mathbf r ^ \prime - \mathbf r | ^ {2} } \delta \left ( \omega - \frac{\mathbf r - \mathbf r ^ \prime }{| \mathbf r - \mathbf r ^ \prime | } \right ) . $$

Here $ \sigma _ {s} ( \mathbf r ) $ is the scattering coefficient (section), $ \sigma ( \mathbf r ) $ is the relaxation coefficient, $ g ( \mu ) $ is the scattering indicatrix, and $ \tau ( \mathbf r ^ \prime , \mathbf r ) $ is the optical length of a path from $ \mathbf r ^ \prime $ to $ \mathbf r $( see [3], [4]). For the construction of estimates with small variances one uses, for example, asymptotic solutions of the adjoint transport equation [4]; the simplest algorithm of this type is the so-called exponential transformation (see [4], [11]). A modification of a local estimate of the flow of particles has been developed (see [3], [4], [11][13], [17], [18]). Using a simulation of a Markov chain (for example, the physical transport process in a certain medium) it is possible to simultaneously obtain dependent estimators of functionals for different values of the parameters; by differentiating the "weight" $ Q _ {n} $ it is sometimes possible to obtain unbiased estimators of the corresponding derivatives (see [4], [12]). This provides an opportunity to use Monte-Carlo methods in solving certain inverse problems [12]. For the solution of problems in transport theory "bifurcation" of the trajectory and analytic averaging is effective [11]. The simulation of trajectories of particles in compound media can sometimes be essentially simplified by the method of a maximal section (see [3][5]).

These are constructed on the basis of the corresponding integral relations. For example, the standard five-point difference approximation for the Laplace equation has the form of the formula for the complete mathematical expectation of a symmetric random walk over a grid with absorption at the boundary (see, for example, [2], [3]). A continuous analogue of this formula is the relation

$$ \tag{2 } u ( P) = \ \frac{\int\limits _ {S ( P) } u ( \mathbf r ( s) ) d s }{\int\limits _ {S ( P) } d s } , $$

where the integral is taken over the surface of a sphere lying entirely within the given domain and with centre at $ P $. Formula (2) and other similar relations provide an opportunity to use isotropic "random walk on spheres" when solving elliptic and parabolic equations (see [24], [4]). Monte-Carlo methods are effective, for example, for estimating the solution of multi-dimensional boundary value problems at a point.

A simulation of Markov branching processes allows one to construct estimates of the solution of certain non-linear equations, for example, the Boltzmann equation in the theory of rarefied gases [3].

References

[1] J. von Neumann, Nat. Bureau Standard Appl. Math. Series , 12 (1951) pp. 36–38
[2] N.P. Buslenko, et al., "The method of statistical trials (the Monte-Carlo method)" , Moscow (1962) (In Russian)
[3] S.M. Ermakov, "Die Monte-Carlo Methode und verwandte Fragen" , Deutsch. Verlag Wissenschaft. (1975) (Translated from Russian)
[4] G.A. Mikhailov, "Some questions in the theory of Monte-Carlo methods" , Novosibirsk (1971) (In Russian)
[5] I.M. Sobol', "Numerical Monte-Carlo methods" , Moscow (1973) (In Russian)
[6] Yu.G. Pollyak, "Probabilistic simulation on computers" , Moscow (1971) (In Russian)
[7] N.S. Bakhvalov, "On optimal convergence estimates for quadrature processes and integration methods of Monte-Carlo type on function classes" , Numerical Methods for Solving Differential and Integral Equations and Quadrature Formulas , Moscow (1964) pp. 5–63 (In Russian)
[8] N.S. Bakhvalov, "An estimate of the remainder term in quadrature formula" Comp. Math. Math. Phys. , 1 : 1 (1961) pp. 68–82 Zh. Vychisl. Mat. i Mat. Fiz. , 1 : 1 (1961) pp. 64–77
[9] N.S. Bakhvalov, "Approximate computation of multiple integrals" Vestnik Moskov. Gos. Univ. Ser. Mat. Mekh. Astronom. Fiz. Khim. , 4 (1959) pp. 3–18 (In Russian)
[10] J.M. Hammersley, D.C. Handscomb, "Monte-Carlo methods" , Methuen (1964)
[11] , Monte-Carlo methods and problems of radiative transport , Moscow (1967)
[12] G.I. Marchuk, et al., "The Monte-Carlo method in atmospheric optics" , Novosibirsk (1976) (In Russian)
[13] J. Spanier, E. Gelbard, "Monte-Carlo principles and neutron transport problems" , Addison-Wesley (1969)
[14] V.V. Chavchanidze, Izv. Akad. Nauk SSSR Ser. Fiz. , 19 : 6 (1955) pp. 629–638
[15] , The penetration of radiation through non-uniformity in protection , Moscow (1968) (In Russian)
[16] A.D. Frank-Kamenetskii, Atomnaya Energiya , 16 : 2 (1964) pp. 119–122
[17] M.H. Kalos, Nuclear Sci. and Eng. , 33 (1968) pp. 284–290
[18] , Monte-Carlo methods and their application. Reports III All-Union Conf. Monte-Carlo Methods , Novosibirsk (1971)
[19] I.M. Gelfand, A.S. Frolov, N.N. Chentsov, "The computation of continuous integrals by the Monte-Carlo method" Izv. Vuz. Ser. Mat. , 5 (1958) pp. 32–45 (In Russian)
[20] A.S. Frolov, N.N. Chentsov, "On the calculation of definite integrals dependent on a parameter by the Monte-Carlo method" USSR Comp. Math. Math. Phys. , 2 : 4 (1962) pp. 802–807 Zh. Vychisl. Mat. i. Mat. Fiz. , 2 : 4 (1962) pp. 714–717
[21] N.M. Korobov, "Number-theoretic methods in applied analysis" , Moscow (1963) (In Russian)
[22] V.S. Vladimirov, "Monte-Carlo methods as applied to the calculation of the lowest eigenvalue and the associated eigenfunction of a linear integral equation" Theor. Probab. Appl. , 1 : 1 (1956) pp. 101–116 Teor. Veroyatnost. i. Primenen. , 1 : 1 (1956) pp. 113–130
[23] J.H. Curtiss, "Monte-Carlo methods for the iteration of linear operators" J. Math. Phys. , 32 : 4 (1954) pp. 209–232
[24] M.E. Muller, "Some continuous Monte-Carlo methods for the Dirichlet problem" Ann. Math. Stat. , 27 : 3 (1956) pp. 569–589

Comments

For an up-to-date account of the use of random processes in solving (numerically) the classical equations of mathematical physics cf. [a2].

References

[a1] P. Bratley, B.L. Fox, L.E. Schrage, "A guide to simulation" , Springer (1987)
[a2] S.M. Ermakov, V.V. Nekrutkin, A.S. Sipin, "Random processes for the classical equations of mathematical physics" , Kluwer (1989) (Translated from Russian)
How to Cite This Entry:
Monte-Carlo method. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Monte-Carlo_method&oldid=47895
This article was adapted from an original article by G.A. Mikhailov (originator), which appeared in Encyclopedia of Mathematics - ISBN 1402006098. See original article