% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Monte Carlo integration}

\section{Introduction}
{\em Monte Carlo integration} is a quadrature (cubature) where the
nodes are chosen randomly. Typically no assumptions are made about the
smoothness of the integrand, not even that it is continuous.

Monte Carlo algorithms are particularly suited for multi-dimensional
integrations where one of the problems is that the integration region,
$\Omega$, often has quite complicated boundary which can not be easily
described by simple functions.  However, it is usually much easier to
find out whether a given point lies within the integration region or not.
Therefore a popular strategy is to create an auxiliary rectangular volume,
$V$, which encompasses the integration volume $\Omega$, and an auxiliary
function which coincides with the integrand inside the volume $\Omega$
and is equal zero outside. Then the integral of the auxiliary function
over the auxiliary volume is equal the original integral.

However, the auxiliary function is generally non-continuous at the
boundary; thus ordinary quadratures which assume continuity of the
integrand are bound to have difficulties here.  One the contrary
the Monte-Carlo quadratures will do just as good (or as bad) as with
continuous integrands.

A typical implementation of a Monte Carlo algorithm integrates the given
function over a rectangular volume, specified by the coordinates of its
"lower-left" and "upper-right" vertices, assuming the user has provided
the encompassing volume with the auxiliary function.

Plain Monte Carlo algorithm distributes points uniformly throughout the
integration region using uncorrelated pseudo-random sequences of points.

Adaptive algorithms, such as VEGAS and MISER, distribute points
non-uniformly in an attempt to reduce integration error using
correspondingly {\em importance} and {\em stratified} sampling.

Yet another strategy to reduce the error is to use correlated quasi-random
sequences.

The GNU Scientific Library, GSL, implements a plain Monte Carlo
integration algorithm; a stratified sampling algorithm, MISER; an
importance sampling algorithm, VEGAS; and a number of quasi-random
generators.

\section{Plain Monte Carlo sampling}
Plain Monte Carlo is a quadrature with random abscissas and equal
weights,
	\begin{equation}
\int_Vf(\ve{x})dV\approx w\sum_{i=1}^Nf(\ve{x}_i)\;,
	\end{equation}
where $\ve{x}$ is a point in the multi-dimensional integration space.
One free parameter, $w$, allows one condition to be satisfied: the
quadrature must integrate exactly a constant function. This gives
$w=V/N$,
	\begin{equation}
\int_Vf(\ve{x})dV \approx \frac{V}{N} \sum_{i=1}^{N}f(\ve{x}_i) \doteq
V\langle f\rangle\;.
	\end{equation}
Under the assumptions of the {\em central limit theorem} the error
of the integration can be estimated as
	\begin{equation}\label{mc:e:ee}
\mathrm{error} = V{\sigma\over \sqrt{N}}\;,
	\end{equation}
where $\sigma$ is the variance of the sample,
	\begin{equation}
\sigma^2 = \langle f^2\rangle-\langle f\rangle^2 \:.
	\end{equation}
The familiar $1/\sqrt{N}$ convergence of a random walk process is quite
slow: to reduce the error by a factor 10 requires 100-fold increase in
the number of sample points.

Expression~(\ref{mc:e:ee}) provides only a statistical estimate of the
error, which is not a strict bound; random sampling may not uncover all
the important features of the function, resulting in an underestimate
of the error.

A simple implementation of the plain Monte Carlo algorithm is shown
in Table~\ref{mc:t:p}.
	\begin{table}[htb]
\caption{Plain Monte Carlo integrator}\label{mc:t:p}
\lstinputlisting{c/monte.carlo/plainmc.c}
	\end{table}


\section{Importance sampling}
Suppose the points are distributed not uniformly but with some density
$\rho(\ve{x})\;$.  That is, the number of points $\Delta n$ in the volume
$\Delta V$ around point $\ve{x}$ is given as
	\begin{equation}
\Delta n = \frac{N}{V}\rho(\ve{x})\Delta V,
	\end{equation}
where $\rho$ is normalised such that $\int_V\rho dV=V$.

The estimate of the integral is then given as
	\begin{equation}
\int_V f(\ve{x}) dV
\approx \sum_{i=1}^{N}f(\ve{x}_i)\Delta V_i
= \sum_{i=1}^{N}f(\ve{x}_i) \frac{V}{N\rho(\ve{x}_i)}
= V\left\langle{f\over\rho}\right\rangle\;,
	\end{equation}
where
	\begin{equation}
\Delta V_i=\frac{V}{N\rho(x_i)}
	\end{equation}
is the volume-per-point at the point $x_i$.

The corresponding variance is now given by
\begin{equation}
\sigma^2 = \left\langle\left({f\over\rho}\right)^2\right\rangle
-\left\langle{f\over\rho}\right\rangle^2\;.
\end{equation}
Apparently if the ratio $f/\rho$ is close to a constant, the variance
is reduced.

It is tempting to take $\rho=\left|f\right|$ and sample directly from the
integrand. However in practice evaluations of the integrand are typically
expensive. Therefore a better strategy is to build an approximate density
in the product form, $\rho(x,y,\dots,z)=\rho_x(x)\rho_y(y)\dots\rho_z(z)$,
and then sample from this approximate density. A popular routine of
this sort is called VEGAS.

%The sampling from a given function can be done using the {\em Metropolis
%algorithm} which shall not be discussed here.

\section{Stratified sampling}

Stratified sampling is a generalisation of the recursive adaptive
integration algorithm to random quadratures in multi-dimensional spaces.

	\begin{table}[htb]
	\caption{Recursive stratified sampling algorithm}
	\label{tb:mc:ss}
\begin{lstlisting}[mathescape=true]
sample $N$ random points with plain Monte Carlo;
estimate the average and the error;
IF the error is acceptable :
	RETURN the average and the error;
ELSE : 
	FOR EACH dimension :
		subdivide the volume in two along the dimension;
		estimate the sub-variances in the two sub-volumes;
	pick the dimension with the largest sub-variance;
	subdivide the volume in two along this dimension;
	dispatch two recursive calls to each of the sub-volumes;
	estimate the grand average and grand error;
	RETURN the grand average and grand error;
\end{lstlisting}
\end{table}

The ordinary ``dividing by two'' strategy does not work for
multi-dimensional integrations as the number of sub-volumes grows way
too fast to keep track of. Instead one estimates along which dimension a
subdivision should bring the most dividends and only subdivides along this
dimension.  Such strategy is called {\em recursive stratified sampling}. A
simple variant of this algorithm is presented in Table~\ref{tb:mc:ss}.

\begin{figure}[htb]
%\centerline{\input{javascript/monte.carlo/strata.tex}}
\centerline{\includegraphics{javascript/monte.carlo/strata.eps}}
\caption{Stratified sample of a discontinuous function\hfill$ $
\centerline{$f(x,y)=\left\{
\mathop{}_{0 \;\mathrm{otherwise}}^{1\;\mathrm{if}\:x^2+y^2<0.9^2}
\right.$} built with the algorithm in Table~\ref{tb:mc:ss}.}
	\label{mc:f:strata}
	\end{figure}

In a stratified sample the points are concentrated in the regions
where the variance of the function is largest, see an illustration in
Figure~\ref{mc:f:strata}.


\section{Quasi-random (low-discrepancy) sampling}

\begin{figure}[tb]
\centerline{
\includegraphics{c/monte.carlo/pseudo.ps}
\includegraphics{c/monte.carlo/quasi.ps}
\includegraphics{c/monte.carlo/halton.ps}
}
\caption{Typical distributions of pseudo-random points (left),
and quasi-random low-discrepancy points: lattice (center) and base-2/3
Halton (right) sequences. The first thousand points are plotted in each
case.}
	\label{fig:pseudo}
	\end{figure}

Pseudo-random sampling has high discrepancy\footnote{discrepancy is a
measure of how unevenly the points are distributed over the region.}: it
typically creates regions with hight density of points and other regions
with low density of points, see an illustration on Figure~\ref{fig:pseudo}
(left).  With pseudo-random sampling there is a finite probability that
all the $N$ points would fall into one half of the region and none into
the other half.

{\em Quasi-random} sequences avoid this phenomenon by distributing
points in a highly correlated manner with a specific requirement of
low discrepancy, see Figure~\ref{fig:pseudo} for an example. Quasi-random
sampling is like a computation on a grid where the grid constant must
not be known in advance as the grid is ever gradually refined and the
points are always distributed uniformly over the region. The computation
can be stopped at any time.

By placing points more evenly than at random, the quasi-random sequences
try to improve on the $1/\sqrt{N}$ convergence rate of pseudo-random
sampling.

The central limit theorem does not apply in this case as the points
are not statistically independent. Therefore the variance can not be
used as an estimate of the error. The error estimation is actually not
trivial. In practice one can employ two different sequences and use the
difference in the resulting integrals as an error estimate.

%Quasi-random sequences can be roughly divided into {\em lattice rules}
%and {\em digital nets}.
% (see e.g. arXiv:1003.4785 [math.NA] and references
%therein).

\subsection{Van der Corput and Halton sequences}
A {\em van der Corput sequence} is a low-discrepancy sequence over the
unit interval. It is constructed by reversing the base-$b$ representation
of the sequence of natural numbers $(1, 2, 3, \dots)$.	For example,
the decimal van der Corput sequence begins as
	\begin{equation}\label{}
0.1, 0.2, 0.3,\dots,0.8, 0.9, 0.01, 0.11, 0.21, 0.31,\dots,
0.91,0.02,0.12,\dots \;.
	\end{equation}

In a base-$b$ representation a natural number $n$ with $s$
digits $\{d_i\mid\:i=1 \ldots s,\; 0\le d_i<b\}$ is given as
	\begin{equation}\label{}
n = \sum_{k=1}^{s} d_k b^{k-1} \;.
	\end{equation}
The corresponding base-$b$ van der Corput number $q_b(n)$ is then given as
	\begin{equation}\label{mc:e:c}
q_b(n) = \sum_{k=1}^{s} d_k b^{-k} \;.
	\end{equation}
Here is a C implementation of this algorithm,
	\begin{lstlisting}
double corput(int n, int base){
   double q=0, bk=(double)1/base;
   while(n>0){ q += (n % base)*bk; n /= base; bk /= base; }
   return q; }
	\end{lstlisting}

The van der Corput numbers of any base are uniformly distributed over
the unit interval.  They also form a dense set in the unit interval:
there exists a subsequence of the van der Corput sequence which converges
to any given real number in $[0,1]$.

The Halton sequence is a generalization of the van der Corput
sequence to $d$-dimensional spaces.  One chooses a set of coprime
bases $b_1,\dots,b_d$ and then for each dimension $i$ generates a
van der Corput sequence with its own base $b_i$.  The $n$-th Halton
$d$-dimentional point $\ve{x}$ in the unit volume is then given as
	\begin{equation}\label{mc:eq:h}
\ve{x}_{b_1,\dots,b_d}(n) = \{q_{b_1}(n),\dots,q_{b_d}(n)\} \;.
	\end{equation}
Here is a C implementation which calls the {\tt corput} function listed
above,
	\begin{lstlisting}
#include<assert.h>
void halton(int n, int d, double *x){
  int base[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67};
  int maxd=19; assert(d <= maxd);
  for(int i=0;i<d;i++) x[i]=corput(n,base[i]); }
	\end{lstlisting}


\subsection{Lattice rules}

In the simplest incarnation a lattice rule can be defined as
follows.  Let the {\em generating vector} $\ve{z}=\{\alpha_i \:|\:
i=1,\dots,d\}$---where $d$ is the dimension of the integration space---be
a set of cleverly chosen irrational numbers.  Then the $n$-th point
(in the unit volume) of the sequence is given as
	\begin{equation}
\ve{x}(n)= \mathrm{frac}(n\ve{z}) \equiv
\{\mathrm{frac}(n\alpha_1),\dots,\mathrm{frac}(n\alpha_d)\}\;,
	\end{equation}
where $\mathrm{frac}(x)$ is the fractional part of $x$.

An implementation of this algorithm in C is given in Table~\ref{mc:t:l}
and an illustration of such sequence is shown on Figure~\ref{fig:pseudo}
(center).

	\begin{table}[htb]
	\caption{Lattice low-disrepancy quasi-random sequence in C.}
	\label{mc:t:l}
	\lstinputlisting{c/monte.carlo/lattice.c} \end{table}

\section{Implementations}

	\begin{table}[htb]
	\caption{Javascript implementation of the stratified sampling algorithm.}
	\label{mc:t:jss}
	\lstset{frame=single,basicstyle=\footnotesize,tabsize=2}
	\lstinputlisting{javascript/monte.carlo/strata.js} \end{table}
