382 lines
16 KiB
Markdown
382 lines
16 KiB
Markdown
# Exercise 5
|
|
|
|
The following integral is to be evaluated comparing different Monte Carlo
|
|
techniques.
|
|
\vspace{20pt}
|
|
$$
|
|
I = \int_0^1 \!\!\! dx \, e^x
|
|
$$
|
|
\vspace{-70pt}
|
|
\begin{figure}
|
|
\hypertarget{fig:exp}{%
|
|
\centering
|
|
\begin{tikzpicture}
|
|
% Integral
|
|
\filldraw [cyclamen!15!white, domain=0:5, variable=\x]
|
|
(0,0) -- plot({\x},{exp(\x/5)}) -- (5,0) -- cycle;
|
|
\draw [cyclamen] (5,0) -- (5,2.7182818);
|
|
\node [below] at (5,0) {1};
|
|
% Axis
|
|
\draw [thick, <-] (0,4) -- (0,0);
|
|
\draw [thick, ->] (-2,0) -- (7,0);
|
|
\node [below right] at (7,0) {$x$};
|
|
\node [above left] at (0,4) {$e^{x}$};
|
|
% Plot
|
|
\draw [domain=-2:7, smooth, variable=\x,
|
|
cyclamen, ultra thick] plot ({\x},{exp(\x/5)});
|
|
\end{tikzpicture}
|
|
\caption{Plot of the integral to be evaluated.}
|
|
}
|
|
\end{figure}
|
|
whose exact value is 1.7182818285...
|
|
|
|
The three most popular Monte Carlo (MC) methods where applied: plain MC, Miser
|
|
and Vegas. Besides being commonly used, these were chosen for also being
|
|
implemented in the GSL libraries `gsl_monte_plain`, `gsl_monte_miser` and
|
|
`gsl_monte_vegas`.
|
|
|
|
|
|
## Plain Monte Carlo
|
|
|
|
When an integral $I$ over a $n-$dimensional space $\Omega$ of volume $V$ of a
|
|
function $f$ has to be evaluated, that is:
|
|
$$
|
|
I = \int\limits_{\Omega} dx \, f(x)
|
|
\with V = \int\limits_{\Omega} dx
|
|
$$
|
|
the simplest MC method approach is to sample $N$ points $x_i$ in $V$ and
|
|
approximate $I$ as:
|
|
$$
|
|
I \approx I_N = \frac{V}{N} \sum_{i=1}^N f(x_i) = V \cdot \avg{f}
|
|
$$
|
|
If $x_i$ are uniformly distributed $I_N \rightarrow I$ for $N \rightarrow +
|
|
\infty$ by the law of large numbers, whereas the integral variance can be
|
|
estimated as:
|
|
$$
|
|
\sigma^2_f = \frac{1}{N - 1}
|
|
\sum_{i = 1}^N \left( f(x_i) - \avg{f} \right)^2
|
|
\et
|
|
\sigma^2_I = \frac{V^2}{N^2} \sum_{i = 1}^N
|
|
\sigma^2_f = \frac{V^2}{N} \sigma^2_f
|
|
$$
|
|
Thus, the error decreases as $1/\sqrt{N}$.
|
|
Unlike in deterministic methods, the error estimate is not a strict bound:
|
|
random sampling may not cover all the important features of the integrand and
|
|
this can result in an underestimation of the error.
|
|
|
|
In this case $f(x) = e^{x}$ and $\Omega = [0,1]$, hence $V = 1$.
|
|
|
|
{#fig:plain-mc-iter}
|
|
|
|
Since the distance from $I$ of $I_N$ is related to $N$, the accuracy of the
|
|
method lies in how many points are generated, namely how many function calls
|
|
are executed when the iterative method is implemented. In @fig:plain-mc-iter
|
|
and @fig:miser-iter, results obtained with the plain MC method are shown in
|
|
red. In @tbl:plain-mc-res, some of them are listed: the estimated integrals
|
|
$I^{\text{oss}}$ are compared to the expected value $I$ and the differences
|
|
between them are given.
|
|
|
|
---------------------------------------------------------------------------
|
|
calls $I^{\text{oss}}$ $\sigma$ diff
|
|
------------------ ------------------ ------------------ ------------------
|
|
\num{5e5} 1.7166435813 0.0006955691 0.0016382472
|
|
|
|
\num{5e6} 1.7181231109 0.0002200309 0.0001587176
|
|
|
|
\num{5e7} 1.7183387184 0.0000695809 0.0000568899
|
|
---------------------------------------------------------------------------
|
|
|
|
Table: Some MC results with three different numbers of function calls.
|
|
Differences between computed and exact values are given in
|
|
diff. {#tbl:plain-mc-res}
|
|
|
|
As can be seen, $\sigma$ is always of the same order of magnitude of diff,
|
|
except for very low numbers of function calls. Even with \num{5e7} calls,
|
|
$I^{\text{oss}}$ still differs from $I$ at the fifth decimal place, meaning
|
|
that this method shows a really slow convergence.
|
|
|
|
The $\sigma$ dependence on the number $C$ of function calls was checked with a
|
|
least square minimization by modeling the data with the function:
|
|
$$
|
|
\sigma = \frac{a}{x^b}
|
|
$$
|
|
The obtained result (@fig:err_fit) confirms the expected result
|
|
$b^{\text{exp}} = 0.5$, with an observed value of $\sim 0.499$.
|
|
Given this dependence, for an error of $10^{-n}$, a number $\propto 10^{2n}$ of
|
|
function calls is needed. To compute an integral within double precision, an
|
|
impossibly large number of $\sigma \sim 10^{32}$ calls is needed, which makes
|
|
this method unpractical for high-precision applications.
|
|
|
|
{#fig:err_fit}
|
|
|
|
|
|
## Stratified sampling
|
|
|
|
In statistics, stratified sampling is a method of sampling from a population
|
|
partitioned into subpopulations. Stratification, indeed, is the process of
|
|
dividing the primary sample into subgroups (strata) before sampling
|
|
within each stratum.
|
|
Given a sample $\{x_j\}_i$ of the $i$-th strata, its mean $\bar{x}_i$ and
|
|
variance ${\sigma^2_x}_i$, are given by:
|
|
$$
|
|
\bar{x}_i = \frac{1}{n_i} \sum_j x_j
|
|
$$
|
|
and
|
|
$$
|
|
\sigma_i^2 = \frac{1}{n_i - 1} \sum_j \left( x_j - \bar{x}_i \right)^2
|
|
\thus
|
|
{\sigma^2_x}_i = \frac{1}{n_i^2} \sum_j \sigma_i^2 = \frac{\sigma_i^2}{n_i}
|
|
$$
|
|
where:
|
|
|
|
- $j$ runs over the points $x_j$ of the sample,
|
|
- $n_i$ is the size of the sample,
|
|
- $\sigma_i^2$ is the variance associated to every point of the $i$-th
|
|
stratum.
|
|
|
|
An estimation of the mean $\bar{x}$ and variance $\sigma_x^2$ for the whole
|
|
population are then given by the stratified sampling as follows:
|
|
$$
|
|
\bar{x} = \frac{1}{N} \sum_i N_i \bar{x}_i \et
|
|
\sigma_x^2 = \sum_i \left( \frac{N_i}{N} \right)^2 {\sigma_x}^2_i
|
|
= \sum_i \left( \frac{N_i}{N} \right)^2 \frac{\sigma^2_i}{n_i}
|
|
$$
|
|
where:
|
|
|
|
- $i$ runs over the strata,
|
|
- $N_i$ is the weight of the $i$-th stratum
|
|
- $N$ is the sum of all strata weights.
|
|
|
|
In practical terms, it can produce a weighted mean that has less variability
|
|
than the arithmetic mean of a simple random sample of the whole population. In
|
|
fact, if measurements within strata have lower standard deviation, the final
|
|
result will have a smaller error in estimation with respect to the one otherwise
|
|
obtained with simple sampling.
|
|
For this reason, stratified sampling is used as a method of variance reduction
|
|
when MC methods are used to estimate population statistics from a known
|
|
population. For examples, see [@ridder17].
|
|
|
|
|
|
### MISER
|
|
|
|
The MISER technique aims at reducing the integration error through the use of
|
|
recursive stratified sampling.
|
|
As stated before, according to the law of large numbers, for a great number of
|
|
extracted points, the estimation of the integral $I$ can be computed as:
|
|
$$
|
|
I= V \cdot \avg{f}
|
|
$$
|
|
Since $V$ is known (in this case, $V = 1$), it is sufficient to estimate
|
|
$\avg{f}$.
|
|
|
|
Consider two disjoint regions $a$ and $b$, such that $a \cup b = \Omega$, in
|
|
which $n_a$ and $n_b$ points are respectively uniformly sampled. Given the
|
|
Monte Carlo estimates of the means $\avg{f}_a$ and $\avg{f}_b$ of those points
|
|
and their variances $\sigma_a^2$ and $\sigma_b^2$, if the weights $N_a$ and
|
|
$N_b$ of $\avg{f}_a$ and $\avg{f}_b$ are chosen unitary, then the variance
|
|
$\sigma^2$ of the combined estimate $\avg{f}$:
|
|
$$
|
|
\avg{f} = \frac{1}{2} \left( \avg{f}_a + \avg{f}_b \right)
|
|
$$
|
|
is given by:
|
|
$$
|
|
\sigma^2 = \frac{\sigma_a^2}{4n_a} + \frac{\sigma_b^2}{4n_b}
|
|
$$
|
|
It can be shown that this variance is minimized by distributing the points such
|
|
that:
|
|
$$
|
|
\frac{n_a}{n_a + n_b} = \frac{\sigma_a}{\sigma_a + \sigma_b}
|
|
$$
|
|
Hence, the smallest error estimate is obtained by allocating sample points in
|
|
proportion to the standard deviation of the function in each sub-region.
|
|
The whole integral estimate and its variance are therefore given by:
|
|
$$
|
|
I = V \cdot \avg{f} \et \sigma_I^2 = V^2 \cdot \sigma^2
|
|
$$
|
|
When implemented, MISER is in fact a recursive method. First, all the possible
|
|
bisections of $\Omega$ are tested and the one which minimizes the combined
|
|
variance of the two sub-regions is selected. In order to speed up the
|
|
algorithm, the variance in the sub-regions is estimated with a fraction of the
|
|
total number of available points (function calls), in GSL it is default set to
|
|
0.1. The remaining points are allocated to the sub-regions using the formula for
|
|
$n_a$ and $n_b$, once the variances are computed.
|
|
|
|
This procedure is then repeated recursively for each of the two half-regions
|
|
from the best bisection. When the allocated calls for a region running out
|
|
(less than 512 in GSL), the method falls back to a plain Monte Carlo.
|
|
The final individual values and their error estimates are then combined upwards
|
|
to give an overall result and an estimate of its error [@sayah19].
|
|
|
|
{#fig:miser-iter}
|
|
|
|
Results for this particular sample are shown in black in @fig:miser-iter and
|
|
compared with the plain MC results (in red). Some of them are listed in
|
|
@tbl:miser-res. Except for the first very little number of calls, the
|
|
improvement with respect to the Plain MC technique is appreciable.
|
|
|
|
The convergence is much faster than a plain MC: at 500'000 function calls, the
|
|
estimate agrees with the exact integral to the fifth decimal place. Once again,
|
|
the standard deviation and the difference share the same magnitude.
|
|
|
|
--------------------------------------------------------------------
|
|
calls $I^{\text{oss}}$ $\sigma$ diff
|
|
----------- ------------------ ------------------ ------------------
|
|
\num{5e5} 1.7182850738 0.0000021829 0.0000032453
|
|
|
|
\num{5e6} 1.7182819143 0.0000001024 0.0000000858
|
|
|
|
\num{5e7} 1.7182818221 0.0000000049 0.0000000064
|
|
--------------------------------------------------------------------
|
|
|
|
Table: MISER results with different numbers of function calls. Differences
|
|
between computed and exact values are given in diff. {#tbl:miser-res}
|
|
|
|
|
|
## Importance sampling
|
|
|
|
In Monte Carlo methods, importance sampling is a technique which samples points
|
|
from distribution whose shape is close to the integrand $f$ itself. This way,
|
|
the points cluster in the regions that make the largest contribution to the
|
|
integral $\int f(x)dx$ and consequently decrease the variance.
|
|
|
|
In a plain MC the points are sampled uniformly, so their probability
|
|
density is given by:
|
|
$$
|
|
g(x) = \frac{1}{V} \quad \forall x \in \Omega
|
|
$$
|
|
and the integral can be written as:
|
|
$$
|
|
I = \int_\Omega dx f(x) = V \int_\Omega f(x) \frac{1}{V}dx
|
|
\approx V \avg{f}
|
|
$$
|
|
More generally, consider a distribution $h(x)$ and similarly do:
|
|
$$
|
|
I
|
|
= \int_\Omega dx f(x)
|
|
= \int_\Omega dx \, \frac{f(x)}{h(x)} \, h(x)
|
|
= \Exp \left[ \frac{f}{h}, h \right]
|
|
$$
|
|
where $\Exp[X, h]$ is the expected value of $X$ wrt $h$. Also note that $h$ has
|
|
to vanish outside $\Omega$ for this to hold.
|
|
To reduce the variance, as anticipated, $h$ must be close to $f$. Assuming they
|
|
are proportional, $h(x) = \alpha |f(x)|$, it follows that:
|
|
$$
|
|
\Exp \left[ \frac{f}{h}, h \right] = \frac{1}{\alpha}
|
|
\et
|
|
\Var \left[ \frac{f}{h}, h \right] = 0
|
|
$$
|
|
For the expected value to give the original $I$, the proportionality constant
|
|
must be taken to be $I^{-1}$, meaning:
|
|
$$
|
|
h(z) = \frac{1}{I}\, |f(z)|
|
|
$$
|
|
The sampling from this $h$ would produce a perfect result with zero variance.
|
|
Of course, this is nonsense: if $I$ is known in advance, there would be no need
|
|
to do a Monte Carlo integration to begin with. Nonetheless, this example serves
|
|
to prove how variance reduction is achieved by sampling from an approximation
|
|
of the integrand.
|
|
|
|
In conclusion, since certain values of $x$ have more impact on $\Exp[f/h, h]$
|
|
than others, these "important" values must be emphasized by sampling them more
|
|
frequently. As a consequence, the estimator variance will be reduced.
|
|
|
|
|
|
### VEGAS
|
|
|
|
The VEGAS algorithm [@lepage78] of G. P. Lepage is based on importance
|
|
sampling. As stated before, it is in practice impossible to sample points from
|
|
the best distribution $h(x)$: only a good approximation can be achieved. The
|
|
VEGAS algorithm attempts this by building a histogram of the function $f$ in
|
|
different subregions with an iterative method, namely:
|
|
|
|
- a fixed number of points (function calls) is generated uniformly in the
|
|
whole region;
|
|
|
|
- the volume $V$ is divided into $N$ intervals of width $\Delta x_i =
|
|
\Delta x \, \forall \, i$, where $N$ is limited by the computer storage
|
|
space available and must be held constant from iteration to iteration.
|
|
(In GSL this defaults to $N = 50$);
|
|
|
|
- each interval is then divided into $m_i + 1$ subintervals, where:
|
|
$$
|
|
m_i = K \frac{\bar{f}_i \Delta x_i}{\sum_j \bar{f}_j \Delta x_j}
|
|
$$
|
|
|
|
where $j$ runs over all the intervals and $\bar{f}_i$ is the average value
|
|
of $f$ in the interval. Hence, $m_i$ is therefore a measure of the
|
|
"importance" of the interval with respect to the others: the higher
|
|
$\bar{f}_i$, the higher $m_i$. The constant $K$, called stiffness, defaults
|
|
to 1.5 in GSL;
|
|
|
|
- as it is desirable to restore the number of intervals to its original value
|
|
$N$, groups of the new intervals must be merged into larger intervals, the
|
|
number of subintervals in each group being constant. The net effect is
|
|
to alter the intervals sizes, while keeping the total number constant, so
|
|
that the smallest intervals occur where $f$ is largest;
|
|
|
|
- the function is integrated with a plain MC method in each interval
|
|
and the sum of the integrals is taken as the $j$-th estimate of $I$. Its
|
|
error is given by the sum of the variances in each interval;
|
|
|
|
- the new grid is further refined in subsequent iterations. By default, the
|
|
number of iterations is 5 in GSL.
|
|
|
|
The final estimate of the integral $I$ and its error
|
|
$\sigma_I$ are made based on weighted average:
|
|
$$
|
|
\avg{I} = \sigma_I^2 \sum_i \frac{I_i}{\sigma_i^2}
|
|
\with
|
|
\frac{1}{\sigma_I^2} = \sum_i \frac{1}{\sigma_i^2}
|
|
$$
|
|
where $I_i$ and $\sigma_i$ are are the integral and standard deviation
|
|
estimated in each iteration.
|
|
|
|
The reliability of the result is asserted by a chi-squared per degree of
|
|
freedom $\chi_r^2$, which should be close to 1 for a good estimation. At a
|
|
given iteration $i$, the $\chi^2_i$ is computed as follows:
|
|
$$
|
|
\chi^2_i = \sum_{j \le i}
|
|
\frac{(I_j - \avg{I})^2}{\sigma_j^2}
|
|
$$
|
|
While performing the iterations, if the value of $\chi_r^2$ exceeds 1.5, the
|
|
routine stops since it is not making progress.
|
|
|
|
Clearly, a better estimation is achieved with a greater number of function
|
|
calls. For this particular sample, the most accurate results are shown in
|
|
@fig:vegas-iter and some of them are listed in @tbl:vegas-res.
|
|
|
|
{#fig:vegas-iter}
|
|
|
|
As can be appreciated in @fig:vegas-iter, the VEGAS algorithm manages to compute the
|
|
integral value more accurately compared to MISER. The $\chi_r^2$ turns out to
|
|
be close enough to 1 to guarantee a good estimation of $I$, goodness which is
|
|
also confirmed by the very small difference shown in @tbl:vegas-res.
|
|
In fact, with a number of \num{5e7} function calls, the difference is
|
|
smaller than \num{1e-10}.
|
|
|
|
----------------------------------------------------------------------------------------------
|
|
calls $I^{\text{oss}}$ $\sigma$ diff $\chi_r^2$
|
|
------------------ ------------------ ------------------ ------------------ ------------------
|
|
\num{5e5} 1.7182818281 0.0000000012 0.0000000004 1.457
|
|
|
|
\num{5e6} 1.7182818284 0.0000000000 0.0000000001 0.632
|
|
|
|
\num{5e7} 1.7182818285 0.0000000000 0.0000000000 0.884
|
|
----------------------------------------------------------------------------------------------
|
|
|
|
Table: Best VEGAS results with different numbers of
|
|
function calls. {#tbl:vegas-res}
|
|
|
|
In conclusion, between a plain Monte Carlo technique, stratified sampling and
|
|
importance sampling, the last turned out to be the most powerful mean to
|
|
obtain a good estimation of the integrand.
|