diff --git a/.gitignore b/.gitignore new file mode 100755 index 000000000..a757647ff --- /dev/null +++ b/.gitignore @@ -0,0 +1,54 @@ +# compiled files: +*.o +*.so +*.a +# temporary files: +*.pyc +*.bak +*.swp +*~ +.*~ +*.old +tmp* +temp* +.#* +\#* +# tex files: +*.log +*.dvi +*.aux +*.blg +*.idx +*.nav +*.out +*.toc +*.snm +*.vrb +# eclipse files: +*.cproject +*.project +# misc: +.DS_Store +Trash +# doconce: +.*_html_file_collection + +*.copyright + +*.copyright + +*.exerinfo + +doc/Projects/2016/Project2/.out.copyright + +doc/Projects/2016/Project2/.project2_2016.copyright + +*.exerinfo + +*.dlog + + + +*.user + +doc/Programs/ParallelizationMPI/dill10 diff --git a/LICENSE b/LICENSE new file mode 100644 index 000000000..2c4afabdb --- /dev/null +++ b/LICENSE @@ -0,0 +1,117 @@ +CC0 1.0 Universal + +Statement of Purpose + +The laws of most jurisdictions throughout the world automatically confer +exclusive Copyright and Related Rights (defined below) upon the creator and +subsequent owner(s) (each and all, an "owner") of an original work of +authorship and/or a database (each, a "Work"). + +Certain owners wish to permanently relinquish those rights to a Work for the +purpose of contributing to a commons of creative, cultural and scientific +works ("Commons") that the public can reliably and without fear of later +claims of infringement build upon, modify, incorporate in other works, reuse +and redistribute as freely as possible in any form whatsoever and for any +purposes, including without limitation commercial purposes. These owners may +contribute to the Commons to promote the ideal of a free culture and the +further production of creative, cultural and scientific works, or to gain +reputation or greater distribution for their Work in part through the use and +efforts of others. + +For these and/or other purposes and motivations, and without any expectation +of additional consideration or compensation, the person associating CC0 with a +Work (the "Affirmer"), to the extent that he or she is an owner of Copyright +and Related Rights in the Work, voluntarily elects to apply CC0 to the Work +and publicly distribute the Work under its terms, with knowledge of his or her +Copyright and Related Rights in the Work and the meaning and intended legal +effect of CC0 on those rights. + +1. Copyright and Related Rights. A Work made available under CC0 may be +protected by copyright and related or neighboring rights ("Copyright and +Related Rights"). Copyright and Related Rights include, but are not limited +to, the following: + + i. the right to reproduce, adapt, distribute, perform, display, communicate, + and translate a Work; + + ii. moral rights retained by the original author(s) and/or performer(s); + + iii. publicity and privacy rights pertaining to a person's image or likeness + depicted in a Work; + + iv. rights protecting against unfair competition in regards to a Work, + subject to the limitations in paragraph 4(a), below; + + v. rights protecting the extraction, dissemination, use and reuse of data in + a Work; + + vi. database rights (such as those arising under Directive 96/9/EC of the + European Parliament and of the Council of 11 March 1996 on the legal + protection of databases, and under any national implementation thereof, + including any amended or successor version of such directive); and + + vii. other similar, equivalent or corresponding rights throughout the world + based on applicable law or treaty, and any national implementations thereof. + +2. Waiver. To the greatest extent permitted by, but not in contravention of, +applicable law, Affirmer hereby overtly, fully, permanently, irrevocably and +unconditionally waives, abandons, and surrenders all of Affirmer's Copyright +and Related Rights and associated claims and causes of action, whether now +known or unknown (including existing as well as future claims and causes of +action), in the Work (i) in all territories worldwide, (ii) for the maximum +duration provided by applicable law or treaty (including future time +extensions), (iii) in any current or future medium and for any number of +copies, and (iv) for any purpose whatsoever, including without limitation +commercial, advertising or promotional purposes (the "Waiver"). Affirmer makes +the Waiver for the benefit of each member of the public at large and to the +detriment of Affirmer's heirs and successors, fully intending that such Waiver +shall not be subject to revocation, rescission, cancellation, termination, or +any other legal or equitable action to disrupt the quiet enjoyment of the Work +by the public as contemplated by Affirmer's express Statement of Purpose. + +3. Public License Fallback. Should any part of the Waiver for any reason be +judged legally invalid or ineffective under applicable law, then the Waiver +shall be preserved to the maximum extent permitted taking into account +Affirmer's express Statement of Purpose. In addition, to the extent the Waiver +is so judged Affirmer hereby grants to each affected person a royalty-free, +non transferable, non sublicensable, non exclusive, irrevocable and +unconditional license to exercise Affirmer's Copyright and Related Rights in +the Work (i) in all territories worldwide, (ii) for the maximum duration +provided by applicable law or treaty (including future time extensions), (iii) +in any current or future medium and for any number of copies, and (iv) for any +purpose whatsoever, including without limitation commercial, advertising or +promotional purposes (the "License"). The License shall be deemed effective as +of the date CC0 was applied by Affirmer to the Work. Should any part of the +License for any reason be judged legally invalid or ineffective under +applicable law, such partial invalidity or ineffectiveness shall not +invalidate the remainder of the License, and in such case Affirmer hereby +affirms that he or she will not (i) exercise any of his or her remaining +Copyright and Related Rights in the Work or (ii) assert any associated claims +and causes of action with respect to the Work, in either case contrary to +Affirmer's express Statement of Purpose. + +4. Limitations and Disclaimers. + + a. No trademark or patent rights held by Affirmer are waived, abandoned, + surrendered, licensed or otherwise affected by this document. + + b. Affirmer offers the Work as-is and makes no representations or warranties + of any kind concerning the Work, express, implied, statutory or otherwise, + including without limitation warranties of title, merchantability, fitness + for a particular purpose, non infringement, or the absence of latent or + other defects, accuracy, or the present or absence of errors, whether or not + discoverable, all to the greatest extent permissible under applicable law. + + c. Affirmer disclaims responsibility for clearing rights of other persons + that may apply to the Work or any use thereof, including without limitation + any person's Copyright and Related Rights in the Work. Further, Affirmer + disclaims responsibility for obtaining any necessary consents, permissions + or other rights required for any use of the Work. + + d. Affirmer understands and acknowledges that Creative Commons is not a + party to this document and has no duty or obligation with respect to this + CC0 or use of the Work. + +For more information, please see + + diff --git a/doc/src/Statistics/Statistics.do.txt b/doc/src/Statistics/Statistics.do.txt new file mode 100644 index 000000000..d7891df5e --- /dev/null +++ b/doc/src/Statistics/Statistics.do.txt @@ -0,0 +1,3146 @@ +TITLE: Computational Physics Lectures: Introduction to Monte Carlo methods +AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University +DATE: today + +!split +===== Monte Carlo methods, plan for the lectures ===== +!bblock +o Intro, MC integration and probability distribution functions (PDFs) +o More on integration, PDFs, MC integration and random walks. +o Random walks and statistical physics. +o Statistical physics and the Ising and Potts models +o Quantum Monte Carlo + +!eblock + + +!split +===== Monte Carlo: Enhances algorithmic thinking! ===== +!bblock +* Be able to generate random variables following a given probability distribution function PDF +* Find a probability distribution function (PDF) +* Sampling rule for accepting a move +* Compute standard deviation and other expectation values +* Techniques for improving errors +!eblock + +!split +===== Domains and probabilities ===== +!bblock +Consider the following simple example, namely the tossing of a dice, resulting in the following possible values +!bt +\begin{equation*} +\{2,3,4,5,6,7,8,9,10,11,12\}. +\end{equation*} +!et +These values are called the *domain*. +To this domain we have the corresponding *probabilities* +!bt +\begin{equation*} +\{1/36,2/36/3/36,4/36,5/36,6/36,5/36,4/36,3/36,2/36,1/36\}. +\end{equation*} +!et +!eblock + +!split +===== Monte Carlo methods, tossing a dice ===== +!bblock +The numbers in the domain are the outcomes of the physical process tossing the dice. +We cannot tell beforehand whether the outcome is 3 or 5 or any other number in this domain. +This defines the randomness of the outcome, or unexpectedness or any other synonimous word which +encompasses the uncertitude of the final outcome. + +The only thing we can tell beforehand +is that say the outcome 2 has a certain probability. +If our favorite hobby is to spend an hour every evening throwing dice and +registering the sequence of outcomes, we will note that the numbers in the above domain +!bt +\begin{equation*} +\{2,3,4,5,6,7,8,9,10,11,12\}, +\end{equation*} +!et +appear in a random order. After 11 throws the results may look like + +!bt +\begin{equation*} +\{10,8,6,3,6,9,11,8,12,4,5\}. +\end{equation*} +!et +!eblock + +!split +===== Stochastic variables ===== +!bblock + +_Random variables are characterized by a domain which contains all possible values that the random value may take. This domain has a corresponding PDF_. +!eblock + +!split +===== Stochastic variables and the main concepts, the discrete case ===== +!bblock +There are two main concepts associated with a stochastic variable. The +*domain* is the set $\mathbb D = \{x\}$ of all accessible values +the variable can assume, so that $X \in \mathbb D$. An example of a +discrete domain is the set of six different numbers that we may get by +throwing of a dice, $x\in\{1,\,2,\,3,\,4,\,5,\,6\}$. + +The *probability distribution function (PDF)* is a function +$p(x)$ on the domain which, in the discrete case, gives us the +probability or relative frequency with which these values of $X$ +occur +!bt +\begin{equation*} +p(x) = \mathrm{Prob}(X=x). +\end{equation*} +!et +!eblock + + +!split +===== Stochastic variables and the main concepts, the continuous case ===== +!bblock +In the continuous case, the PDF does not directly depict the +actual probability. Instead we define the probability for the +stochastic variable to assume any value on an infinitesimal interval +around $x$ to be $p(x)dx$. The continuous function $p(x)$ then gives us +the *density* of the probability rather than the probability +itself. The probability for a stochastic variable to assume any value +on a non-infinitesimal interval $[a,\,b]$ is then just the integral + +!bt +\begin{equation*} +\mathrm{Prob}(a\leq X\leq b) = \int_a^b p(x)dx. +\end{equation*} +!et +Qualitatively speaking, a stochastic variable represents the values of +numbers chosen as if by chance from some specified PDF so that the +selection of a large set of these numbers reproduces this PDF. +!eblock + +!split +===== The cumulative probability ===== +!bblock +Of interest to us is the *cumulative probability +distribution function* (_CDF_), $P(x)$, which is just the probability +for a stochastic variable $X$ to assume any value less than $x$ +!bt +\begin{equation*} +P(x)=\mathrm{Prob(}X\leq x\mathrm{)} = +\int_{-\infty}^x p(x^{\prime})dx^{\prime}. +\end{equation*} +!et +The relation between a CDF and its corresponding PDF is then + +!bt +\begin{equation*} +p(x) = \frac{d}{dx}P(x). +\end{equation*} +!et +!eblock + +!split +===== Properties of PDFs ===== +!bblock + +There are two properties that all PDFs must satisfy. The first one is +positivity (assuming that the PDF is normalized) + +!bt +\begin{equation*} +0 \leq p(x) \leq 1. +\end{equation*} +!et +Naturally, it would be nonsensical for any of the values of the domain +to occur with a probability greater than $1$ or less than $0$. Also, +the PDF must be normalized. That is, all the probabilities must add up +to unity. The probability of ``anything'' to happen is always unity. For +both discrete and continuous PDFs, this condition is +!bt +\begin{align*} +\sum_{x_i\in\mathbb D} p(x_i) & = 1,\\ +\int_{x\in\mathbb D} p(x)\,dx & = 1. +\end{align*} +!et +!eblock + +!split +===== Important distributions, the uniform distribution ===== +!bblock +The first one +is the most basic PDF; namely the uniform distribution +!bt +\begin{equation} +p(x) = \frac{1}{b-a}\theta(x-a)\theta(b-x), +label{eq:unifromPDF} +\end{equation} +!et +with +!bt +\begin{equation*} +\begin{array}{ll} +\theta(x)=0 & x<0 \\ +\theta(x)=\frac{1}{b-a} & \in [a,b]. +\end{array} +\end{equation*} +!et +The normal distribution with $b=1$ and $a=0$ is used to generate random numbers. +!eblock + +!split +===== Gaussian distribution ===== +!bblock +The second one is the Gaussian Distribution +!bt +\begin{equation*} +p(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp{(-\frac{(x-\mu)^2}{2\sigma^2})}, +\end{equation*} +!et +with mean value $\mu$ and standard deviation $\sigma$. If $\mu=0$ and $\sigma=1$, it is normally called the _standard normal distribution_ +!bt +\begin{equation*} +p(x) = \frac{1}{\sqrt{2\pi}} \exp{(-\frac{x^2}{2})}, +\end{equation*} +!et + +The following simple Python code plots the above distribution for different values of $\mu$ and $\sigma$. +!bc pyscpro +import numpy as np +from math import acos, exp, sqrt +from matplotlib import pyplot as plt +from matplotlib import rc, rcParams +import matplotlib.units as units +import matplotlib.ticker as ticker +rc('text',usetex=True) +rc('font',**{'family':'serif','serif':['Gaussian distribution']}) +font = {'family' : 'serif', + 'color' : 'darkred', + 'weight' : 'normal', + 'size' : 16, + } +pi = acos(-1.0) +mu0 = 0.0 +sigma0 = 1.0 +mu1= 1.0 +sigma1 = 2.0 +mu2 = 2.0 +sigma2 = 4.0 + +x = np.linspace(-20.0, 20.0) +v0 = np.exp(-(x*x-2*x*mu0+mu0*mu0)/(2*sigma0*sigma0))/sqrt(2*pi*sigma0*sigma0) +v1 = np.exp(-(x*x-2*x*mu1+mu1*mu1)/(2*sigma1*sigma1))/sqrt(2*pi*sigma1*sigma1) +v2 = np.exp(-(x*x-2*x*mu2+mu2*mu2)/(2*sigma2*sigma2))/sqrt(2*pi*sigma2*sigma2) +plt.plot(x, v0, 'b-', x, v1, 'r-', x, v2, 'g-') +plt.title(r'{\bf Gaussian distributions}', fontsize=20) +plt.text(-19, 0.3, r'Parameters: $\mu = 0$, $\sigma = 1$', fontdict=font) +plt.text(-19, 0.18, r'Parameters: $\mu = 1$, $\sigma = 2$', fontdict=font) +plt.text(-19, 0.08, r'Parameters: $\mu = 2$, $\sigma = 4$', fontdict=font) +plt.xlabel(r'$x$',fontsize=20) +plt.ylabel(r'$p(x)$ [MeV]',fontsize=20) + +# Tweak spacing to prevent clipping of ylabel +plt.subplots_adjust(left=0.15) +plt.savefig('gaussian.pdf', format='pdf') +plt.show() +!ec +!eblock + + +!split +===== Exponential distribution ===== +!bblock +Another important distribution in science is the exponential distribution +!bt +\begin{equation*} +p(x) = \alpha\exp{-(\alpha x)}. +\end{equation*} +!et +!eblock + +!split +===== Expectation values ===== +!bblock +Let $h(x)$ be an arbitrary continuous function on the domain of the stochastic +variable $X$ whose PDF is $p(x)$. We define the *expectation value* +of $h$ with respect to $p$ as follows + +!bt +\begin{equation} +\langle h \rangle_X \equiv \int\! h(x)p(x)\,dx +label{eq:expectation_value_of_h_wrt_p} +\end{equation} +!et +Whenever the PDF is known implicitly, like in this case, we will drop +the index $X$ for clarity. +A particularly useful class of special expectation values are the +*moments*. The $n$-th moment of the PDF $p$ is defined as +follows +!bt +\begin{equation*} +\langle x^n \rangle \equiv \int\! x^n p(x)\,dx +\end{equation*} +!et +!eblock + +!split +===== Stochastic variables and the main concepts, mean values ===== +!bblock +The zero-th moment $\langle 1\rangle$ is just the normalization condition of +$p$. The first moment, $\langle x\rangle$, is called the *mean* of $p$ +and often denoted by the letter $\mu$ +!bt +\begin{equation*} +\langle x\rangle = \mu \equiv \int x p(x)dx, +\end{equation*} +!et +for a continuous distribution and +!bt +\begin{equation*} +\langle x\rangle = \mu \equiv \frac{1}{N}\sum_{i=1}^N x_i p(x_i), +\end{equation*} +!et +for a discrete distribution. +Qualitatively it represents the centroid or the average value of the +PDF and is therefore simply called the expectation value of $p(x)$. +!eblock + +!split +===== Stochastic variables and the main concepts, central moments, the variance ===== +!bblock + +A special version of the moments is the set of *central moments*, the n-th central moment defined as +!bt +\begin{equation*} +\langle (x-\langle x\rangle )^n\rangle \equiv \int\! (x-\langle x\rangle)^n p(x)\,dx +\end{equation*} +!et +The zero-th and first central moments are both trivial, equal $1$ and +$0$, respectively. But the second central moment, known as the +*variance* of $p$, is of particular interest. For the stochastic +variable $X$, the variance is denoted as $\sigma^2_X$ or $\mathrm{Var}(X)$ +!bt +\begin{align*} +\sigma^2_X &=\mathrm{Var}(X) = \langle (x-\langle x\rangle)^2\rangle = +\int (x-\langle x\rangle)^2 p(x)dx\\ +& = \int\left(x^2 - 2 x \langle x\rangle^{2} +\langle x\rangle^2\right)p(x)dx\\ +& = \langle x^2\rangle\rangle - 2 \langle x\rangle\langle x\rangle + \langle x\rangle^2\\ +& = \langle x^2 \rangle - \langle x\rangle^2 +\end{align*} +!et +The square root of the variance, $\sigma =\sqrt{\langle (x-\langle x\rangle)^2\rangle}$ is called the +_standard deviation_ of $p$. It is the RMS (root-mean-square) +value of the deviation of the PDF from its mean value, interpreted +qualitatively as the ``spread'' of $p$ around its mean. +!eblock + + + +!split +===== First Illustration of the Use of Monte-Carlo Methods, integration ===== +!bblock +With this definition of a random variable and its associated PDF, +we attempt now a clarification of the Monte-Carlo strategy by using the +evaluation of an integral +as our example. + +In discussion on numerical integration we went through standard methods for evaluating an integral like + +!bt +\begin{equation*} + I=\int_0^1 f(x)dx\approx \sum_{i=1}^N\omega_if(x_i), +\end{equation*} +!et +where $\omega_i$ are the weights determined by the specific integration method +(like Simpson's method) with $x_i$ the given mesh points. +To give you a feeling of how we are to evaluate the above integral using Monte-Carlo, +we employ here the crudest possible approach. Later on we will present +slightly more refined approaches. +This crude approach consists in setting all weights equal 1, $\omega_i=1$. That corresponds to the +rectangle method +!bt +\begin{equation*} + I=\int_a^bf(x) dx \approx h\sum_{i=1}^N f(x_{i-1/2}), +\end{equation*} +!et +where $f(x_{i-1/2})$ is the midpoint value of $f$ for a given value $x_{i-1/2}$. +!eblock + +!split +===== First Illustration of the Use of Monte-Carlo Methods, integration ===== +!bblock +Setting $h=(b-a)/N$ where $b=1$, $a=0$, we can then rewrite the above integral as +!bt +\begin{equation*} + I=\int_0^1 f(x)dx\approx \frac{1}{N}\sum_{i=1}^Nf(x_{i-1/2}), +\end{equation*} +!et +where $x_{i-1/2}$ are the midpoint values of $x$. +Introducing the concept of the average of the function $f$ for a given PDF $p(x)$ as +!bt +\begin{equation*} + \langle f \rangle = \frac{1}{N}\sum_{i=1}^Nf(x_i)p(x_i), +\end{equation*} +!et +and identify $p(x)$ with the uniform distribution, viz +$ p(x)=1$ when $x\in [0,1]$ and zero for all other values of $x$. +The integral is +is then the average of $f$ over the interval $x \in [0,1]$ + +!bt +\begin{equation*} + I=\int_0^1 f(x)dx\approx \langle f \rangle. +\end{equation*} +!et +!eblock + +!split +===== First Illustration of the Use of Monte-Carlo Methods, variance in integration ===== +!bblock +In addition to the average value $\langle f \rangle$ the other +important quantity in a +Monte-Carlo calculation is the variance $\sigma^2$ and +the standard deviation $\sigma$. We define first the variance +of the integral with $f$ for a uniform distribution in the interval +$x \in [0,1]$ to be +!bt +\begin{equation*} + \sigma^2_f=\frac{1}{N}\sum_{i=1}^N(f(x_i)-\langle f\rangle)^2p(x_i), +\end{equation*} +!et +and inserting the uniform distribution this yields +!bt +\begin{equation*} + \sigma^2_f=\frac{1}{N}\sum_{i=1}^Nf(x_i)^2- + \left(\frac{1}{N}\sum_{i=1}^Nf(x_i)\right)^2, +\end{equation*} +!et +or +!bt +\begin{equation*} + \sigma^2_f=\left(\langle f^2\rangle - + \langle f \rangle^2\right). +\end{equation*} +!et +!eblock + +!split +===== Monte-Carlo integration, meaning of variance ===== +!bblock +The variance is nothing but a measure of the extent to +which $f$ deviates from its average over the region of integration. +The standard deviation is defined as the square root of the variance. +If we consider the above results for +a fixed value of $N$ as a measurement, +we could recalculate the +above average and variance for a series of different measurements. +If each such measumerent produces a set of averages for the +integral $I$ denoted $\langle f\rangle_l$, we have for $M$ measurements +that the integral is given by +!bt +\begin{equation*} + \langle I \rangle_M=\frac{1}{M}\sum_{l=1}^{M}\langle f\rangle_l. +\end{equation*} +!et +!eblock + +!split +===== First Illustration of the Use of Monte-Carlo Methods, integration ===== +!bblock +If we can consider the probability of +correlated events to be zero, we can rewrite +the variance of these series of measurements as (equating $M=N$) +!bt +\begin{equation} + \sigma^2_N\approx \frac{1}{N}\left(\langle f^2\rangle - + \langle f \rangle^2\right)=\frac{\sigma^2_f}{N}. +label{eq:sigmaN} +\end{equation} +!et +We note that the standard deviation is proportional to the inverse square root of +the number of measurements +!bt +\begin{equation*} + \sigma_N \sim \frac{1}{\sqrt{N}}. +\end{equation*} +!et +!eblock + + +!split +===== Important aspects of Monte-Carlo Methods ===== +!bblock +*The aim of Monte Carlo calculations is to have $\sigma_N$ as small as possible after $N$ samples. * +The results from one sample represents, +since we are using concepts from statistics, +a 'measurement'. +!eblock + +!split +===== Why Monte Carlo integration? ===== +!bblock +The scaling in the previous equation +is clearly unfavorable compared even with the +trapezoidal rule. We saw that the trapezoidal +rule carries a truncation error +!bt +\[ +\mathrm{error}\sim O(h^2), +\] +!et +with $h$ the step length. +In general, methods based on a Taylor expansion such as the trapezoidal +rule or Simpson's rule have a truncation +error which goes like $\sim O(h^k)$, with $k \ge 1$. +Recalling that the step size is defined as $h=(b-a)/N$, we have an +error which goes like +!bt +\[ +\mathrm{error}\sim N^{-k}. +\] +!et +!eblock + +!split +===== Why Monte Carlo integration? ===== +!bblock +Monte Carlo integration is more efficient in higher dimensions. +To see this, let us assume that our integration volume is a hypercube +with side $L$ and dimension $d$. This cube contains hence +$N=(L/h)^d$ points and therefore the error in the result scales as +$N^{-k/d}$ for the traditional methods. + +The error in the Monte carlo integration is +however independent of $d$ and scales as +!bt +\[ +\mathrm{error}\sim 1/\sqrt{N}. +\] +!et +_Always_! + +Comparing this error with that of the traditional methods, shows that +Monte Carlo integration is more efficient than an algorithm with error in powers of $k$ +when +!bt +\[ +d>2k. +\] +!et +!eblock + +!split +===== Why Monte Carlo integration? Example ===== +!bblock +In order to expose this, consider the definition of the quantum mechanical energy +of a system consisting of 10 particles in three dimensions. The energy is the expectation value +of the Hamiltonian $H$ and reads +!bt +\begin{equation*} + E=\frac{\int d\mathbf{R}_1d\mathbf{R}_2\dots d\mathbf{R}_N + \Psi^{\ast}(\mathbf{R_1},\mathbf{R}_2,\dots,\mathbf{R}_N) + H(\mathbf{R_1},\mathbf{R}_2,\dots,\mathbf{R}_N) + \Psi(\mathbf{R_1},\mathbf{R}_2,\dots,\mathbf{R}_N)} + {\int d\mathbf{R}_1d\mathbf{R}_2\dots d\mathbf{R}_N + \Psi^{\ast}(\mathbf{R_1},\mathbf{R}_2,\dots,\mathbf{R}_N) + \Psi(\mathbf{R_1},\mathbf{R}_2,\dots,\mathbf{R}_N)}, + \end{equation*} +!et +where $\Psi$ is the wave function of the system and $\mathbf{R}_i$ are the coordinates +of each particle. If we want to compute the above integral +using for example Gaussian quadrature and use for example ten mesh +points for the ten particles, we need to compute a ten-dimensional integral with a total of $10^{30}$ mesh points. +As an amusing exercise, assume that you have access to today's fastest computer with a theoretical peak +capacity of more than one Petaflops, that is $10^{15}$ floating point operations per second. Assume also that every mesh point +corresponds to one floating operation per second. Estimate then the time needed to compute this integral with +a traditional method like Gaussian quadrature and compare this number with the estimated lifetime of the +universe, $T\approx 4.7 \times 10^{17}$s. Do you have the patience to wait? +!eblock + + +!split +===== Monte Carlo integration, simple example ===== +!bblock +We end this first part with a discussion of +a brute force Monte Carlo program which integrates +!bt +\begin{equation*} + \int_0^1dx\frac{4}{1+x^2} = \pi, +\end{equation*} +!et +where the input is the desired number of Monte Carlo samples. +!eblock + +!split +===== Monte Carlo integration, simple example ===== +!bblock +What we are doing is to employ a random number generator to obtain numbers +$x_i$ in the interval $[0,1]$ through a call to one of the +library functions +$ran0$, $ran1$, $ran2$ or $ran3$ which generate random numbers in the +interval $x\in [0,1]$. +These functions will be discussed in the next section. +Here we simply employ these functions in order +to generate a random variable. +All random number generators produce pseudo-random +numbers in the interval $[0,1]$ using the so-called uniform +probability distribution +$p(x)$ defined as + +!bt +\begin{equation*} + p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x), +\end{equation*} +!et +with $a=0$ og $b=1$ and where $\Theta$ is the standard Heaviside +function or simply the step function. +!eblock + + +!split +===== Monte Carlo integration, simple example ===== +!bblock +If we have a general interval $[a,b]$, we can still +use these random number generators through a change of variables + +!bt +\begin{equation*} + z=a+(b-a)x, +\end{equation*} +!et +with $x$ in the interval $x\in [0,1]$. +!eblock +!split +===== Monte Carlo integration, simple example ===== +!bblock +The present approach to the above integral is often called 'crude' or +'Brute-Force' Monte-Carlo. +Later on in this chapter we will study refinements to this +simple approach. The reason is that a random generator +produces +points that are distributed +in a homogenous way in the interval $[0,1]$. +If our function is peaked around certain values of $x$, +we may end +up sampling function values where +$f(x)$ is small or near zero. Better schemes which reflect the +properties of the function to be integrated are thence needed. +!eblock + +!split +===== Monte Carlo integration, algorithm ===== +!bblock +The algorithm is as follows + + * Choose the number of Monte Carlo samples $N$. + + * Perform a loop over $N$ and for each step generate a random number $x_i$ in the interval $[0,1]$ through a call to a random number generator. + + * Use this number to evaluate $f(x_i)$. + + * Evaluate the contributions to the mean value and the standard deviation for each loop. + + * After $N$ samples calculate the final mean value and the standard deviation. +!eblock + +!split +===== Monte Carlo integration, simple example, the program ===== +!bblock +!bc cppcod +#include +#include +using namespace std; + +// Here we define various functions called by the main program +// this function defines the function to integrate + +double func(double x); + +// Main function begins here +int main() +{ + int n; + double MCint, MCintsqr2, fx, Variance; + cout << "Read in the number of Monte-Carlo samples" << endl; + cin >> n; + MCint = MCintsqr2=0.; + double invers_period = 1./RAND_MAX; // initialise the random number generator + srand(time(NULL)); // This produces the so-called seed in MC jargon +// evaluate the integral with the a crude Monte-Carlo method + for ( int i = 1; i <= n; i++){ + // obtain a floating number x in [0,1] + double x = double(rand())*invers_period; + fx = func(x); + MCint += fx; + MCintsqr2 += fx*fx; + } + MCint = MCint/((double) n ); + MCintsqr2 = MCintsqr2/((double) n ); + double variance=MCintsqr2-MCint*MCint; +// final output + cout << " variance= " << variance << " Integral = " << MCint << " Exact= " << M_PI << endl; +} // end of main program +// this function defines the function to integrate +double func(double x) +{ + double value; + value = 4/(1.+x*x); + return value; +} // end of function to evaluate +!ec +!eblock + + +!split +===== Monte Carlo integration, simple example and the results ===== +!bblock + +|-----------------------------------------| +| $N$ | $I$ | $\sigma_N$ | +|-----r------------l------------l---------| +| 10 | 3.10263E+00 | 3.98802E-01 | +| 100 | 3.02933E+00 | 4.04822E-01 | +| 1000 | 3.13395E+00 | 4.22881E-01 | +| 10000 | 3.14195E+00 | 4.11195E-01 | +| 100000 | 3.14003E+00 | 4.14114E-01 | +| 1000000 | 3.14213E+00 | 4.13838E-01 | +| 10000000 | 3.14177E+00 | 4.13523E-01 | +| $10^{9}$ | 3.14162E+00 | 4.13581E-01 | +|-----------------------------------------| + +We note that as $N$ increases, the integral itself never reaches more than an agreement +to the fourth or fifth digit. The variance also oscillates around its exact value +$4.13581E-01$. +!eblock + + +!split +===== Testing against the trapezoidal rule for a one-dimensional integral ===== +!bblock +The following simple Python code, with pertaining plot shows the relative error for the above integral using a brute force Monte Carlo approach and the trapezoidal rule. Running the python code shows that the trapezoidal rule is clearly superior in this case. With importance sampling and multi-dimensional integrals, the Monte Carl method takes over again. +!bc pyscpro +from matplotlib import pyplot as plt +from math import exp, acos, log10 +import numpy as np +import random + +# function for the trapezoidal rule +def TrapezoidalRule(a,b,f,n): + h = (b-a)/float(n) + s = 0 + x = a + for i in range(1,n,1): + x = x+h + s = s+ f(x) + s = 0.5*(f(a)+f(b))+s + return h*s +# function to perform the Monte Carlo calculations +def MonteCarloIntegration(f,n): + sum = 0 +# Define the seed for the rng + random.seed() + for i in range (1, n, 1): + x = random.random() + sum = sum +f(x) + return sum/n + +# function to compute +def function(x): + return 4/(1+x*x) + +# Integration limits for the Trapezoidal rule +a = 0.0; b = 1.0 +exact = acos(-1.0) +# set up the arrays for plotting the relative error +log10n = np.zeros(6); Trapez = np.zeros(6); MCint = np.zeros(6); +# find the relative error as function of integration points +for i in range(1, 6): + npts = 10**(i+1) + log10n[i] = log10(npts) + Trapez[i] = log10(abs((TrapezoidalRule(a,b,function,npts)-exact)/exact)) + MCint[i] = log10(abs((MonteCarloIntegration(function,npts)-exact)/exact)) +plt.plot(log10n, Trapez ,'b-',log10n, MCint,'g-') +plt.axis([1,6,-14.0, 0.0]) +plt.xlabel('$\log_{10}(n)$') +plt.ylabel('Relative error') +plt.title('Relative errors for Monte Carlo integration and Trapezoidal rule') +plt.legend(['Trapezoidal rule', 'Brute force Monte Carlo integration'], loc='best') +plt.savefig('mcintegration.pdf') +plt.show() + +!ec +!eblock + + +!split +===== Second example, particles in a box ===== +!bblock + +We give here an example of how a system evolves towards +a well defined equilibrium state. + +Consider a box divided into two equal halves separated by a wall. +At the beginning, time $t=0$, there are $N$ particles on the left +side. A small hole in the wall is then opened and one particle +can pass through the hole per unit time. + +After some time the system reaches its equilibrium state with +equally many particles in both halves, $N/2$. +Instead of determining complicated initial conditions for a system +of $N$ particles, we model the system by a simple statistical model. +In order to simulate this system, which may consist of $N \gg 1$ particles, +we assume that all particles in the left half have equal probabilities +of going to the right half. +We introduce the label $n_l$ to denote the + number of particles at every time on the left side, and $n_r=N-n_l$ for those +on the right side. +!eblock + + +!split +===== Second example, particles in a box ===== +!bblock +The probability for a move to the right during a time step $\Delta t$ +is $n_l/N$. The algorithm for simulating this problem may then look +like this + + * Choose the number of particles $N$. + + b* Make a loop over time, where the maximum time (or maximum number of steps) should be larger than the number of particles $N$. + + * For every time step $\Delta t$ there is a probability $n_l/N$ for a move to the right. Compare this probability with a random number $x$. + + * If $ x \le n_l/N$, decrease the number of particles in the left half by one, i.e., $n_l=n_l-1$. Else, move a particle from the right half to the left, i.e., $n_l=n_l+1$. + + * Increase the time by one unit (the external loop). +!eblock + +!split +===== Second example, particles in a box ===== +!bblock +In this case, a Monte Carlo sample corresponds to one time unit +$\Delta t$. + +The following simple C/C++-program illustrates this model. +!bc cppcod +// Particles in a box +#include +#include +#include +#include "lib.h" +using namespace std; + +ofstream ofile; +int main(int argc, char* argv[]) +{ + char *outfilename; + int initial_n_particles, max_time, time, random_n, nleft; + long idum; + // Read in output file, abort if there are too few command-line arguments + if( argc <= 1 ){ + cout << "Bad Usage: " << argv[0] << + " read also output file on same line" << endl; + exit(1); + } + else{ + outfilename=argv[1]; + } + ofile.open(outfilename); + // Read in data + cout << "Initial number of particles = " << endl ; + cin >> initial_n_particles; + // setup of initial conditions + nleft = initial_n_particles; + max_time = 10*initial_n_particles; + idum = -1; + // sampling over number of particles + for( time=0; time <= max_time; time++){ + random_n = ((int) initial_n_particles*ran0(&idum)); + if ( random_n <= nleft){ + nleft -= 1; + } + else{ + nleft += 1; + } + ofile << setiosflags(ios::showpoint | ios::uppercase); + ofile << setw(15) << time; + ofile << setw(15) << nleft << endl; + } + return 0; +} // end main function +!ec + +!eblock + + +!split +===== Second example, particles in a box, discussion ===== +!bblock +If we denote $\langle n_l \rangle $ as the number of particles in the left +half as a time average after *equilibrium is reached*, +we can define the standard +deviation as + +!bt +\begin{equation} + \sigma =\sqrt{\langle n_l^2 \rangle-\langle n_l \rangle^2}. +\end{equation} +!et + +This problem has also an analytic solution to which we can compare +our numerical simulation. +!eblock + +!split +===== Second example, particles in a box, discussion ===== +!bblock +If $n_l(t)$ is the number of particles in the left half after +$t$ moves, the change in $n_l(t)$ in the time interval $\Delta t$ +is + +!bt +\begin{equation*} + \Delta n=\left(\frac{N-n_l(t)}{N}-\frac{n_l(t)}{N}\right)\Delta t, +\end{equation*} +!et +and assuming that $n_l$ and $t$ are continuous variables we arrive at +!bt +\begin{equation*} +\frac{dn_l(t)}{dt}=1-\frac{2n_l(t)}{N}, +\end{equation*} +!et +whose solution is +!bt +\begin{equation*} + n_l(t)=\frac{N}{2}\left(1+e^{-2t/N}\right), +\end{equation*} +!et +with the initial condition $n_l(t=0)=N$. Note that we have assumed $n$ to be a continuous variable. Obviously, particles are discrete objects. +!eblock + +!split +===== Simple demonstration using python ===== +!bblock +The following simple Python code implements the above algorithm for particles in a box and plots the final number of particles in each part of the box. +!bc pyscpro +#!/usr/bin/env python +from matplotlib import pyplot as plt +from math import exp +import numpy as np +import random + +# initial number of particles +N0 = 1000 +MaxTime = 10*N0 +values = np.zeros(MaxTime) +time = np.zeros(MaxTime) +random.seed() +# initial number of particles in left half +nleft = N0 +for t in range (0, MaxTime, 1): + if N0*random.random() <= nleft: + nleft -= 1 + else: + nleft += 1 + time[t] = t + values[t] = nleft + +# Finally we plot the results +plt.plot(time, values,'b-') +plt.axis([0,MaxTime, N0/4, N0]) +plt.xlabel('$t$') +plt.ylabel('$N$') +plt.title('Number of particles in left half') +plt.savefig('box.pdf') +plt.show() +!ec +The produced figure shows the development of this system as +function of time steps. We note that for $N=1000$ +after roughly $2000$ time steps, +the system has reached the equilibrium state. There are however noteworthy +fluctuations around equilibrium. +!eblock + +!split +===== Brief Summary ===== + +In essence the Monte Carlo method contains the following ingredients + + * A PDF which characterizes the system + + * Random numbers which are generated so as to cover in an as uniform as possible way on the unity interval [0,1]. + + * A sampling rule + + * An error estimation + + * Techniques for improving the errors + +!split +===== Probability Distribution Functions ===== +!bblock + +The following table collects properties of probability distribution functions. +In our notation we reserve the label $p(x)$ for the probability of a certain event, +while $P(x)$ is the cumulative probability. + + +|--------------------------------------------------------------------------------------------------------------------------------------| +| | Discrete PDF | Continuous PDF | +|---------------------l-------------------------------------------c-------------------------------------------c------------------------| +| Domain | $\left\{x_1, x_2, x_3, \dots, x_N\right\}$ | $[a,b]$ | +| Probability | $p(x_i)$ | $p(x)dx$ | +| Cumulative | $P_i=\sum_{l=1}^ip(x_l)$ | $P(x)=\int_a^xp(t)dt$ | +| Positivity | $ 0\le p(x_i)\le 1$ | $ p(x) \ge 0$ | +| Positivity | $ 0\le P_i\le 1$ | $ 0\le P(x)\le 1$ | +| Monotonic | $P_i\ge P_j$ if $x_i\ge x_j$ | $P(x_i)\ge P(x_j)$ if $x_i\ge x_j$ | +| Normalization | $P_N=1$ | $P(b)=1$ | +|--------------------------------------------------------------------------------------------------------------------------------------| + +!eblock + + +!split +===== Probability Distribution Functions ===== +!bblock +With a PDF we can compute expectation values of selected quantities such as + +!bt +\begin{equation*} + \langle x^k\rangle=\frac{1}{N}\sum_{i=1}^{N}x_i^kp(x_i), +\end{equation*} +!et +if we have a discrete PDF or + +!bt +\begin{equation*} + \langle x^k\rangle=\int_a^b x^kp(x)dx, +\end{equation*} +!et +in the case of a continuous PDF. We have already defined the mean value $\mu$ +and the variance $\sigma^2$. +!eblock + +!split +===== The three famous Probability Distribution Functions ===== +!bblock + +There are at least three PDFs which one may encounter. These are the + +_Uniform distribution_ +!bt +\begin{equation*} +p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x), +\end{equation*} +!et +yielding probabilities different from zero in the interval $[a,b]$. + +_The exponential distribution_ +!bt +\begin{equation*} +p(x)=\alpha \exp{(-\alpha x)}, +\end{equation*} +!et +yielding probabilities different from zero in the interval $[0,\infty)$ and with mean value +!bt +\begin{equation*} +\mu = \int_0^{\infty}xp(x)dx=\int_0^{\infty}x\alpha \exp{(-\alpha x)}dx=\frac{1}{\alpha}, +\end{equation*} +!et +!eblock +with variance +!bt +\begin{equation*} +\sigma^2=\int_0^{\infty}x^2p(x)dx-\mu^2 = \frac{1}{\alpha^2}. +\end{equation*} +!et + +!split +===== Probability Distribution Functions, the normal distribution ===== +!bblock +Finally, we have the so-called univariate normal distribution, or just the _normal distribution_ +!bt +\begin{equation*} +p(x)=\frac{1}{b\sqrt{2\pi}}\exp{\left(-\frac{(x-a)^2}{2b^2}\right)} +\end{equation*} +!et +with probabilities different from zero in the interval $(-\infty,\infty)$. +The integral $\int_{-\infty}^{\infty}\exp{\left(-(x^2\right)}dx$ appears in many calculations, its value +is $\sqrt{\pi}$, a result we will need when we compute the mean value and the variance. +The mean value is +!bt +\begin{equation*} + \mu = \int_0^{\infty}xp(x)dx=\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}x \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx, +\end{equation*} +!et +which becomes with a suitable change of variables +!bt +\begin{equation*} + \mu =\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(a+b\sqrt{2}y)\exp{-y^2}dy=a. +\end{equation*} +!et +!eblock + +!split +===== Probability Distribution Functions, the normal distribution ===== +!bblock +Similarly, the variance becomes +!bt +\begin{equation*} + \sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}(x-\mu)^2 \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx, +\end{equation*} +!et +and inserting the mean value and performing a variable change we obtain + +!bt +\begin{equation*} + \sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(b\sqrt{2}y)^2\exp{\left(-y^2\right)}dy= +\frac{2b^2}{\sqrt{\pi}}\int_{-\infty}^{\infty}y^2\exp{\left(-y^2\right)}dy, +\end{equation*} +!et +and performing a final integration by parts we obtain the well-known result $\sigma^2=b^2$. +It is useful to introduce the standard normal distribution as well, defined by $\mu=a=0$, viz. a distribution +centered around zero and with a variance $\sigma^2=1$, leading to + +!bt +\begin{equation} + p(x)=\frac{1}{\sqrt{2\pi}}\exp{\left(-\frac{x^2}{2}\right)}. +\end{equation} +!et +!eblock + +!split +===== Probability Distribution Functions, the cumulative distribution ===== +!bblock + +The exponential and uniform distributions have simple cumulative functions, +whereas the normal distribution does not, being proportional to the so-called +error function $erf(x)$, given by + +!bt +\begin{equation*} +P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x\exp{\left(-\frac{t^2}{2}\right)}dt, +\end{equation*} +!et +which is difficult to evaluate in a quick way. +!eblock + + +!split +===== Probability Distribution Functions, other important distribution ===== +!bblock + +Some other PDFs which one encounters often in the natural sciences are the binomial distribution +!bt +\begin{equation*} + p(x) = \left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} \hspace{0.5cm}x=0,1,\dots,n, +\end{equation*} +!et +where $y$ is the probability for a specific event, such as the tossing of a coin or moving left or right +in case of a random walker. Note that $x$ is a discrete stochastic variable. + +The sequence of binomial trials is characterized by the following definitions + + * Every experiment is thought to consist of $N$ independent trials. + + * In every independent trial one registers if a specific situation happens or not, such as the jump to the left or right of a random walker. + + * The probability for every outcome in a single trial has the same value, for example the outcome of tossing (either heads or tails) a coin is always $1/2$. +!eblock + +!split +===== Probability Distribution Functions, the binomial distribution ===== +!bblock + +In order to compute the mean and variance we need to recall Newton's binomial +formula +!bt +\begin{equation*} + (a+b)^m=\sum_{n=0}^m \left(\begin{array}{c} m \\ n\end{array}\right)a^nb^{m-n}, +\end{equation*} +!et +which can be used to show that + +!bt +\begin{equation*} +\sum_{x=0}^n\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} = (y+1-y)^n = 1, +\end{equation*} +!et +the PDF is normalized to one. +The mean value is +!bt +\begin{equation*} +\mu = \sum_{x=0}^n x\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} = +\sum_{x=0}^n x\frac{n!}{x!(n-x)!}y^x(1-y)^{n-x}, +\end{equation*} +!et +resulting in +!bt +\begin{equation*} +\mu = +\sum_{x=0}^n x\frac{(n-1)!}{(x-1)!(n-1-(x-1))!}y^{x-1}(1-y)^{n-1-(x-1)}, +\end{equation*} +!et +which we rewrite as + +!bt +\begin{equation*} +\mu=ny\sum_{\nu=0}^n\left(\begin{array}{c} n-1 \\ \nu\end{array}\right)y^{\nu}(1-y)^{n-1-\nu} =ny(y+1-y)^{n-1}=ny. +\end{equation*} +!et +!eblock +The variance is slightly trickier to get. It reads $\sigma^2=ny(1-y)$. + +!split +===== Probability Distribution Functions, Poisson's distribution ===== +!bblock + +Another important distribution with discrete stochastic variables $x$ is +the Poisson model, which resembles the exponential distribution and reads +!bt +\begin{equation*} + p(x) = \frac{\lambda^x}{x!} e^{-\lambda} \hspace{0.5cm}x=0,1,\dots,;\lambda > 0. +\end{equation*} +!et +In this case both the mean value and the variance are easier to calculate, + +!bt +\begin{equation*} +\mu = \sum_{x=0}^{\infty} x \frac{\lambda^x}{x!} e^{-\lambda} = \lambda e^{-\lambda}\sum_{x=1}^{\infty} +\frac{\lambda^{x-1}}{(x-1)!}=\lambda, +\end{equation*} +!et +and the variance is $\sigma^2=\lambda$. +!eblock + + + +!split +===== Probability Distribution Functions, Poisson's distribution ===== +!bblock +An example of applications of the Poisson distribution could be the counting +of the number of $\alpha$-particles emitted from a radioactive source in a given time interval. +In the limit of $n\rightarrow \infty$ and for small probabilities $y$, the binomial distribution +approaches the Poisson distribution. Setting $\lambda = ny$, with $y$ the probability for an event in +the binomial distribution we can show that + +!bt +\begin{equation*} +\lim_{n\rightarrow \infty}\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} e^{-\lambda}=\sum_{x=1}^{\infty}\frac{\lambda^x}{x!} e^{-\lambda}. +\end{equation*} +!et +!eblock + +!split +===== Meet the covariance! ===== +!bblock +An important quantity in a statistical analysis is the so-called covariance. + +Consider the set $\{X_i\}$ of $n$ +stochastic variables (not necessarily uncorrelated) with the +multivariate PDF $P(x_1,\dots,x_n)$. The *covariance* of two +of the stochastic variables, $X_i$ and $X_j$, is defined as follows + +!bt +\begin{align} +\mathrm{Cov}(X_i,\,X_j) & = \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle \\ +&=\int\cdots\int (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)P(x_1,\dots,x_n)\,dx_1\dots dx_n, +label{eq:def_covariance} +\end{align} +!et +with +!bt +\begin{equation*} +\langle x_i\rangle = +\int\cdots\int x_i P(x_1,\dots,x_n)\,dx_1\dots dx_n. +\end{equation*} +!et +!eblock + + + +!split +===== Meet the covariance in matrix disguise ===== +!bblock +If we consider the above covariance as a matrix +!bt +\[ +C_{ij} =\mathrm{Cov}(X_i,\,X_j), +\] +!et +then the diagonal elements are just the familiar +variances, $C_{ii} = \mathrm{Cov}(X_i,\,X_i) = \mathrm{Var}(X_i)$. It turns out that +all the off-diagonal elements are zero if the stochastic variables are +uncorrelated. +!eblock + +!split +===== Meet the covariance, uncorrelated events ===== +!bblock + +This is easy to show, keeping in mind the linearity of +the expectation value. Consider the stochastic variables $X_i$ and +$X_j$, ($i\neq j$) +!bt +\begin{align*} +\mathrm{Cov}(X_i,\,X_j) &= \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle\\ +&=\langle x_i x_j - x_i\langle x_j\rangle - \langle x_i\rangle x_j + \langle x_i\rangle\langle x_j\rangle\rangle\\ +&=\langle x_i x_j\rangle - \langle x_i\langle x_j\rangle\rangle - \langle \langle x_i\rangle x_j \rangle + +\langle \langle x_i\rangle\langle x_j\rangle\rangle\\ +&=\langle x_i x_j\rangle - \langle x_i\rangle\langle x_j\rangle - \langle x_i\rangle\langle x_j\rangle + +\langle x_i\rangle\langle x_j\rangle\\ +&=\langle x_i x_j\rangle - \langle x_i\rangle\langle x_j\rangle +\end{align*} +!et +If $X_i$ and $X_j$ are independent, we get +!bt +\[ +\langle x_i x_j\rangle = +\langle x_i\rangle\langle x_j\rangle=\mathrm{Cov}(X_i, X_j) = 0\ \ (i\neq j). +\] +!et +!eblock + + +!split +===== Numerical experiments and the covariance ===== +!bblock + +Now that we have constructed an idealized mathematical framework, let +us try to apply it to empirical observations. Examples of relevant +physical phenomena may be spontaneous decays of nuclei, or a purely +mathematical set of numbers produced by some deterministic +mechanism. It is the latter we will deal with, using so-called pseudo-random +number generators. In general our observations will contain only a limited set of +observables. We remind the reader that +a *stochastic process* is a process that produces sequentially a +chain of values +!bt +\begin{equation*} +\{x_1, x_2,\dots\,x_k,\dots\}. +\end{equation*} +!et +!eblock + + +!split +===== Numerical experiments and the covariance ===== +!bblock +We will call these +values our *measurements* and the entire set as our measured +*sample*. The action of measuring all the elements of a sample +we will call a stochastic *experiment* (since, operationally, +they are often associated with results of empirical observation of +some physical or mathematical phenomena; precisely an experiment). We +assume that these values are distributed according to some +PDF $p_X^{\phantom X}(x)$, where $X$ is just the formal symbol for the +stochastic variable whose PDF is $p_X^{\phantom X}(x)$. Instead of +trying to determine the full distribution $p$ we are often only +interested in finding the few lowest moments, like the mean +$\mu_X^{\phantom X}$ and the variance $\sigma_X^{\phantom X}$. +!eblock + + +!split +===== Numerical experiments and the covariance, actual situations ===== +!bblock +In practical situations however, a sample is always of finite size. Let that +size be $n$. The expectation value of a sample $\alpha$, the _sample mean_, is then defined as follows +!bt +\begin{equation*} +\langle x_{\alpha} \rangle \equiv \frac{1}{n}\sum_{k=1}^n x_{\alpha,k}. +\end{equation*} +!et +The *sample variance* is: +!bt +\begin{equation*} +\mathrm{Var}(x) \equiv \frac{1}{n}\sum_{k=1}^n (x_{\alpha,k} - \langle x_{\alpha} \rangle)^2, +\end{equation*} +!et +with its square root being the *standard deviation of the sample*. +!eblock + + +!split +===== Numerical experiments and the covariance, our observables ===== +!bblock +You can think of the above observables as a set of quantities which define +a given experiment. This experiment is then repeated several times, say $m$ times. +The total average is then +!bt +\begin{equation} +\langle X_m \rangle= \frac{1}{m}\sum_{\alpha=1}^mx_{\alpha}=\frac{1}{mn}\sum_{\alpha, k} x_{\alpha,k}, +label{eq:exptmean} +\end{equation} +!et +where the last sums end at $m$ and $n$. +The total variance is +!bt +\begin{equation*} +\sigma^2_m= \frac{1}{mn^2}\sum_{\alpha=1}^m(\langle x_{\alpha} \rangle-\langle X_m \rangle)^2, +\end{equation*} +!et +which we rewrite as +!bt +\begin{equation} +\sigma^2_m=\frac{1}{m}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle). +label{eq:exptvariance} +\end{equation} +!et +!eblock + +!split +===== Numerical experiments and the covariance, the sample variance ===== +!bblock + +We define also the sample variance $\sigma^2$ of all $mn$ individual experiments as +!bt +\begin{equation} +\sigma^2=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n (x_{\alpha,k}-\langle X_m \rangle)^2. +label{eq:sampleexptvariance} +\end{equation} +!et + + + +These quantities, being known experimental values or the results from our calculations, +may differ, in some cases +significantly, from the similarly named +exact values for the mean value $\mu_X$, the variance $\mathrm{Var}(X)$ +and the covariance $\mathrm{Cov}(X,Y)$. +!eblock + +!split +===== Numerical experiments and the covariance, central limit theorem ===== +!bblock + +The central limit theorem states that the PDF $\tilde{p}(z)$ of +the average of $m$ random values corresponding to a PDF $p(x)$ +is a normal distribution whose mean is the +mean value of the PDF $p(x)$ and whose variance is the variance +of the PDF $p(x)$ divided by $m$, the number of values used to compute $z$. + +The central limit theorem leads then to the well-known expression for the +standard deviation, given by +!bt +\begin{equation*} + \sigma_m= +\frac{\sigma}{\sqrt{m}}. +\end{equation*} +!et + +In many cases the above estimate for the standard deviation, in particular if correlations are strong, may be too simplistic. We need therefore a more precise defintion of the error and the variance in our results. +!eblock + +!split +===== Definition of Correlation Functions and Standard Deviation ===== +!bblock +Our estimate of the true average $\mu_{X}$ is the sample mean $\langle X_m \rangle$ + +!bt +\begin{equation*} +\mu_{X}^{\phantom X} \approx X_m=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n x_{\alpha,k}. +\end{equation*} +!et + + +We can then use Eq. (ref{eq:exptvariance}) +!bt +\begin{equation*} +\sigma^2_m=\frac{1}{mn^2}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle), +\end{equation*} +!et +and rewrite it as +!bt +\begin{equation*} +\sigma^2_m=\frac{\sigma^2}{n}+\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k +#include +#include +#include +using namespace std; +// output file as global variable +ofstream ofile; + +// Main function begins here +int main(int argc, char* argv[]) +{ + int n; + char *outfilename; + + cin >> n; + double MCint = 0.; double MCintsqr2=0.; + double invers_period = 1./RAND_MAX; // initialise the random number generator + srand(time(NULL)); // This produces the so-called seed in MC jargon + // Compute the variance and the mean value of the uniform distribution + // Compute also the specific values x for each cycle in order to be able to + // the covariance and the correlation function + // Read in output file, abort if there are too few command-line arguments + if( argc <= 2 ){ + cout << "Bad Usage: " << argv[0] << + " read also output file and number of cycles on same line" << endl; + exit(1); + } + else{ + outfilename=argv[1]; + } + ofile.open(outfilename); + // Get the number of Monte-Carlo samples + n = atoi(argv[2]); + double *X; + X = new double[n]; + for (int i = 0; i < n; i++){ + double x = double(rand())*invers_period; + X[i] = x; + MCint += x; + MCintsqr2 += x*x; + } + double Mean = MCint/((double) n ); + MCintsqr2 = MCintsqr2/((double) n ); + double STDev = sqrt(MCintsqr2-Mean*Mean); + double Variance = MCintsqr2-Mean*Mean; +// Write mean value and standard deviation + cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl; + + // Now we compute the autocorrelation function + double *autocor; autocor = new double[n]; + for (int j = 0; j < n; j++){ + double sum = 0.0; + for (int k = 0; k < (n-j); k++){ + sum += (X[k]-Mean)*(X[k+j]-Mean); + } + autocor[j] = sum/Variance/((double) n ); + ofile << setiosflags(ios::showpoint | ios::uppercase); + ofile << setw(15) << setprecision(8) << j; + ofile << setw(15) << setprecision(8) << autocor[j] << endl; + } + ofile.close(); // close output file + return 0; +} // end of main program +!ec +!eblock + + +!split +===== Correlation function and which random number generators should I use ===== +!bblock +The following Python code plots the results for the correlation function from the above program. +!bc pyscpro +import numpy as np +from matplotlib import pyplot as plt +# Load in data file +data = np.loadtxt("datafiles/autocor.dat") +# Make arrays containing x-axis and binding energies as function of A +x = data[:,0] +corr = data[:,1] +plt.plot(x, corr ,'ro') +plt.axis([0,1000,-0.2, 1.1]) +plt.xlabel(r'$d$') +plt.ylabel(r'$C_d$') +plt.title(r'autocorrelation function for RNG') +plt.savefig('autocorr.pdf') +plt.show() + +!ec +!eblock + + + +!split +======= Which RNG should I use? ======= +!bblock +* In the library files lib.cpp and lib.h we have included four popular RNGs taken from the widely used textbook "Numerical Recipes":"http://numerical.recipes/". These are called ran0, ran1, ran2 and ran3. +* C++ has a class called _random_. The "random class":"http://www.cplusplus.com/reference/random/" contains a large selection of RNGs and is highly recommended. Some of these RNGs have very large periods making it thereby very safe to use these RNGs in case one is performing large calculations. In particular, the "Mersenne twister random number engine":"http://www.cplusplus.com/reference/random/mersenne_twister_engine/" has a period of $2^{19937}$. + +!eblock + + +!split +===== How to use the Mersenne generator ===== +!bblock +The following part of a c++ code (from project 4) sets up the uniform distribution for $x\in [0,1]$. +!bc cppcod +/* + +// You need this +#include + +// Initialize the seed and call the Mersienne algo +std::random_device rd; +std::mt19937_64 gen(rd()); +// Set up the uniform distribution for x \in [[0, 1] +std::uniform_real_distribution RandomNumberGenerator(0.0,1.0); + +// Now use the RNG +int ix = (int) (RandomNumberGenerator(gen)*NSpins); +!ec +!eblock + + + +!split +======= Improved Monte Carlo Integration ======= +!bblock +We have presented a simple brute force approach +to integration with the Monte Carlo method. There we sampled +over a given number of points distributed uniformly in the interval +$[0,1]$ +!bt +\begin{equation*} + I=\int_0^1 f(x)dx=\langle f\rangle. +\end{equation*} +!et +Here we introduce two important steps which in most cases improve +upon the above simple brute force approach with the uniform distribution, namely +* change of variables and +* importance sampling +!eblock + +!split +===== Change of Variables ===== +!bblock +The starting point is always the uniform distribution + +!bt +\begin{equation*} +p(x)dx=\left\{\begin{array}{cc} dx & 0 \le x \le 1\\ + 0 & else\end{array}\right. +\end{equation*} +!et +with $p(x)=1$ and +satisfying + +!bt +\begin{equation*} + \int_{-\infty}^{\infty}p(x)dx=1. +\end{equation*} +!et +All random number generators use the uniform distribution to generate numbers $x\in [0,1]$. +!eblock + +!split +===== Change of Variables ===== +!bblock +When we attempt a +transformation to a new variable +$x\rightarrow y$ +we have to conserve the probability +!bt +\begin{equation*} + p(y)dy=p(x)dx, +\end{equation*} +!et +which for the uniform distribution implies + +!bt +\begin{equation*} + p(y)dy=dx. +\end{equation*} +!et +Let us assume that $p(y)$ is a PDF different from the uniform +PDF $p(x)=1$ with $x \in [0,1]$. +If we integrate the last expression we arrive at + +!bt +\begin{equation*} + x(y)=\int_0^y p(y')dy', +\end{equation*} +!et +which is nothing but the cumulative distribution of $p(y)$, i.e., +!bt +\begin{equation*} + x(y)=P(y)=\int_0^y p(y')dy'. +\end{equation*} +!et +!eblock + + +!split +===== Transformed Uniform Distribution ===== +!bblock +Suppose we have the general uniform distribution + +!bt +\begin{equation*} +p(y)dy=\left\{\begin{array}{cc} \frac{dy}{b-a} & a \le y \le b\\ + 0 & else\end{array}\right. +\end{equation*} +!et +If we wish to relate this distribution to the one in the interval +$x \in [0,1]$ +we have +!bt +\begin{equation*} + p(y)dy=\frac{dy}{b-a}=dx, +\end{equation*} +!et +and integrating we obtain the cumulative function + +!bt +\begin{equation*} + x(y)=\int_a^y \frac{dy'}{b-a}, +\end{equation*} +!et +yielding + +!bt +\begin{equation*} + y=a+(b-a)x, +\end{equation*} +!et +a well-known result! +!eblock + + +!split +===== Exponential Distribution ===== +!bblock +Assume that + +!bt +\begin{equation*} + p(y)=\exp{(-y)}, +\end{equation*} +!et +which is the exponential distribution, important for the analysis +of e.g., radioactive decay. Again, +$p(x)$ is given by the uniform distribution with +$x \in [0,1]$, and +with the assumption that the probability is conserved we have + +!bt +\begin{equation*} + p(y)dy=\exp{(-y)}dy=dx, +\end{equation*} +!et +which yields after integration + +!bt +\begin{equation*} + x(y)=P(y)=\int_0^y \exp{(-y')}dy'=1-\exp{(-y)}, +\end{equation*} +!et +or + +!bt +\begin{equation*} + y(x)=-\ln{(1-x)}. +\end{equation*} +!et +!eblock + +!split +===== Exponential Distribution ===== +!bblock +This gives us the new random variable $y$ in the domain +$y \in [0,\infty)$ +determined through the random variable $x \in [0,1]$ generated by +functions like $ran0$. + +This means that if we can factor out +$\exp{(-y)}$ from an integrand we may have + +!bt +\begin{equation*} + I=\int_0^{\infty}F(y)dy=\int_0^{\infty}\exp{(-y)}G(y)dy +\end{equation*} +!et +which we rewrite as + +!bt +\begin{equation*} + \int_0^{\infty}\exp{(-y)}G(y)dy= + \int_0^{1}G(y(x))dx\approx + \frac{1}{N}\sum_{i=1}^NG(y(x_i)), +\end{equation*} +!et +where $x_i$ is a random number in the interval $[0,1]$. +!eblock + +!split +===== Exponential Distribution ===== +!bblock +We have changed the integration limits in the second integral, since we have performed a change of variables. Since we have used the uniform distribution defined for $x\in [0,1]$, the integration limits change to $0$ and $1$. The variable $y$ is now a function of $x$. +Note also that in practical implementations, our random number generators for the +uniform distribution never return exactly 0 or 1, but we may come very close. + +The algorithm for the last example is rather simple. +In the function which sets up the integral, we simply need +to call one of the random number generators +like $ran0$, $ran1$, $ran2$ or $ran3$ in order to obtain numbers +in the interval [0,1]. We obtain $y$ by the taking the logarithm of +$(1-x)$. Our calling function which sets up the new random +variable $y$ may then include statements like +!bc cppcod +..... +idum=-1; +x=ran0(&idum); +y=-log(1.-x); +..... +!ec +!eblock + +!split +===== Normal Distribution ===== +!bblock +For the normal distribution, expressed here as +!bt +\begin{equation*} + g(x,y)=\exp{(-(x^2+y^2)/2)}dxdy. +\end{equation*} +!et +it is rather difficult to find an inverse since the cumulative +distribution is given by the error function $erf(x)$ + +!bt +\begin{equation*} + \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}}\int_{0}^x e^{-t^2} dt. +\end{equation*} +!et +Both "c++":"http://www.cplusplus.com/reference/cmath/erfc/" and "Fortran":"https://gcc.gnu.org/onlinedocs/gfortran/ERFC.html" have this function as intrinsic ones. +!eblock + +!split +===== Normal Distribution ===== +!bblock +We obviously would like to avoid computing an integral everytime we need a random variable. +If we however switch to polar coordinates, we have +for $x$ and $y$ + +!bt +\begin{equation*} + r=\left(x^2+y^2\right)^{1/2} \hspace{1cm} + \theta =tan^{-1}\frac{x}{y}, +\end{equation*} +!et +resulting in + +!bt +\begin{equation*} + g(r,\theta)=r\exp{(-r^2/2)}drd\theta, +\end{equation*} +!et +where the angle $\theta$ could be given by a uniform +distribution in the region $[0,2\pi]$. +Following example 1 above, this implies simply +multiplying random numbers +$x\in [0,1]$ by $2\pi$. +!eblock + +!split +===== Normal Distribution ===== +!bblock +The variable +$r$, defined for $r \in [0,\infty)$ needs to be related to +to random numbers $x'\in [0,1]$. To achieve that, we introduce a new variable + +!bt +\begin{equation*} + u=\frac{1}{2}r^2, +\end{equation*} +!et +and define a PDF +!bt +\begin{equation*} + \exp{(-u)}du, +\end{equation*} +!et +with $u\in [0,\infty)$. +Using the results from example 2 for the exponential distribution, we have + +!bt +\begin{equation*} + u=-\ln{(1-x')}, +\end{equation*} +!et +where $x'$ is a random number generated for $x'\in [0,1]$. +!eblock + +!split +===== Normal Distribution ===== +!bblock +With +!bt +\begin{equation*} + x=r\cos{(\theta)}=\sqrt{2u}\cos{(\theta)}, +\end{equation*} +!et +and + +!bt +\begin{equation*} + y=r\sin{(\theta)}=\sqrt{2u}\sin{(\theta)}, +\end{equation*} +!et +we can obtain new random numbers $x,y$ through +!bt +\begin{equation*} + x=\sqrt{-2\ln{(1-x')}}\cos{(\theta)}, +\end{equation*} +!et +and +!bt +\begin{equation*} + y=\sqrt{-2\ln{(1-x')}}\sin{(\theta)}, +\end{equation*} +!et +with $x'\in [0,1]$ and $\theta \in 2\pi [0,1]$. +!eblock + +!split +===== Normal Distribution ===== +!bblock +A function which yields such random numbers for the normal +distribution would include statements like +!bc cppcod +..... +idum=-1; +radius=sqrt(-2*ln(1.-ran0(idum))); +theta=2*pi*ran0(idum); +x=radius*cos(theta); +y=radius*sin(theta); +..... +!ec +!eblock + + +!split +===== Importance Sampling ===== +!bblock +With the aid of the above variable transformations we address now +one of the most widely used approaches to Monte Carlo integration, +namely importance sampling. + +Let us assume that $p(y)$ is a PDF whose behavior resembles that of a function +$F$ defined in a certain interval $[a,b]$. The normalization condition is + +!bt +\begin{equation*} + \int_a^bp(y)dy=1. +\end{equation*} +!et +We can rewrite our integral as + +!bt +\begin{equation*} + I=\int_a^b F(y) dy =\int_a^b p(y)\frac{F(y)}{p(y)} dy. + label{eq:impsampl1} +\end{equation*} +!et +!eblock + +!split +===== Importance Sampling ===== +!bblock +Since random numbers are generated for the uniform distribution $p(x)$ +with $x\in [0,1]$, we need to perform a change of variables $x\rightarrow y$ +through +!bt +\begin{equation*} + x(y)=\int_a^y p(y')dy', +\end{equation*} +!et +where we used +!bt +\begin{equation*} + p(x)dx=dx=p(y)dy. +\end{equation*} +!et +If we can invert $x(y)$, we find $y(x)$ as well. +!eblock + + +!split +===== Importance Sampling ===== +!bblock +With this change of variables we can express the integral of +Eq. (ref{eq:impsampl1}) as +!bt +\begin{equation*} + I=\int_a^b p(y)\frac{F(y)}{p(y)} dy=\int_{\tilde{a}}^{\tilde{b}}\frac{F(y(x))}{p(y(x))} dx, +\end{equation*} +!et +meaning that a Monte Carlo evaluation of the above integral gives +!bt +\begin{equation*} +\int_{\tilde{a}}^{\tilde{b}}\frac{F(y(x))}{p(y(x))} dx= +\frac{1}{N}\sum_{i=1}^N\frac{F(y(x_i))}{p(y(x_i))}. +\end{equation*} +!et +!eblock + +!split +===== Importance Sampling ===== +!bblock +Note the change in integration limits from $a$ and $b$ to $\tilde{a}$ and +$\tilde{b}$. +The advantage of such a change of variables in case $p(y)$ follows +closely $F$ is that the integrand becomes smooth and we can sample +over relevant values for the integrand. It is however not trivial +to find such a function $p$. +The conditions on $p$ which allow us to perform these transformations +are + + * $p$ is normalizable and positive definite, + + * it is analytically integrable and + + * the integral is invertible, allowing us thereby to express a new variable in terms of the old one. +!eblock + +!split +===== Importance Sampling ===== +!bblock +The variance is now with the definition +!bt +\begin{equation*} + \tilde{F}=\frac{F(y(x))}{p(y(x))}, +\end{equation*} +!et +given by +!bt +\begin{equation*} + \sigma^2=\frac{1}{N}\sum_{i=1}^N\left(\tilde{F}\right)^2- +\left(\frac{1}{N}\sum_{i=1}^N\tilde{F}\right)^2. +label{eq:standard_is} +\end{equation*} +!et +!eblock + +!split +===== Importance Sampling ===== +!bblock +The algorithm for this procedure is + +* Use the uniform distribution to find the random variable $y$ in the interval [0,1]. The function $p(x)$ is a user provided PDF. + +* Evaluate thereafter +!bt +\begin{equation*} +I=\int_a^b F(x) dx =\int_a^b p(x)\frac{F(x)}{p(x)} dx, +\end{equation*} +!et + by rewriting +!bt +\begin{equation*} + \int_a^b p(x)\frac{F(x)}{p(x)} dx = + \int_{\tilde{a}}^{\tilde{b}}\frac{F(x(y))}{p(x(y))} dy, +\end{equation*} +!et +since + +!bt +\begin{equation*} + \frac{dy}{dx}=p(x). +\end{equation*} +!et + * Perform then a Monte Carlo sampling for +!bt +\begin{equation*} +\int_{\tilde{a}}^{\tilde{b}}\frac{F(x(y))}{p(x(y))} dy\approx \frac{1}{N}\sum_{i=1}^N\frac{F(x(y_i))}{p(x(y_i))}, +\end{equation*} +!et +with $y_i\in [0,1]$, +* and evaluate the variance as well. +!eblock + + +!split +===== Importance Sampling, a simple example ===== +!bblock +Let us look again at the integral +!bt +\[ + I=\int_0^1 F(x) dx = \int_0^1 \frac{1}{1+x^2} dx = \frac{\pi}{4}. +\] +!et +We choose the following PDF (which follows closely the function to +integrate) +!bt +\[ + p(x)=\frac{1}{3}\left(4-2x\right) \hspace{1cm} \int_0^1p(x)dx=1, +\] +!et +resulting in +!bt +\[ + \frac{F(0)}{p(0)}=\frac{F(1)}{p(1)}=\frac{3}{4}. +\] +!et +Check that it fullfils the requirements of a PDF! +We perform then the change of variables (via the Cumulative function) +!bt +\[ + y(x)=\int_0^x p(x')dx'=\frac{1}{3}x\left(4-x\right), +\] +!et +or +!bt +\[ + x=2-\left(4-3y\right)^{1/2} +\] +!et +We have that when $y=0$ then $x=0$ and when $y=1$ we have $x=1$. +!eblock + +!split +===== Importance Sampling, a simple example, a simple plot ===== +!bblock + +!bc pyscpro +from matplotlib import pyplot as plt +from math import exp, acos, log10 +import numpy as np + +# function to integrate +def function(x): + return 1.0/(1+x*x) + +# new probability +def newfunction(x): + return 0.3333333*(4.0-2*x) + +Dim = 100 +x = np.linspace(0.0,1.0,Dim) +f = np.zeros(Dim) +g = np.zeros(Dim) +for i in xrange(Dim): + f[i] = function(x[i]) + g[i] = newfunction(x[i]) + +plt.plot(x, f ,'b-',x, g,'g-') +plt.axis([0,1,0.5, 1.5]) +plt.xlabel('$x$') +plt.ylabel('Functions') +plt.title('Similarities between functions') +plt.legend(['Integrand', 'New integrand'], loc='best') +plt.savefig('newprobability.pdf') +plt.show() + +!ec +!eblock + + + +!split +===== Importance Sampling, a simple example, the code part ===== +!bblock +!bc cppcod +// evaluate the integral with importance sampling + for ( int i = 1; i <= n; i++){ + x = ran0(&idum); // random numbers in [0,1] + y = 2 - sqrt(4-3*x); // new random numbers + fy=3*func(y)/(4-2*y); // weighted function + int_mc += fy; + sum_sigma += fy*fy; + } + int_mc = int_mc/((double) n ); + sum_sigma = sum_sigma/((double) n ); + variance=(sum_sigma-int_mc*int_mc); +!ec +!eblock + + + +!split +===== Importance Sampling, a simple example, and the results ===== +!bblock +The suffix $cr$ stands for the brute force approach +while $is$ stands for the use of importance sampling. +All calculations use ran0 as function to generate the uniform distribution. + +|-----------------------------------------------------------------------| +|$N$ | $I_{cr}$ | $\sigma_{cr}$ | $I_{is}$ | $\sigma_{is}$ | +|-----------------------------------------------------------------------| +| 10000 | 3.13395E+00 | 4.22881E-01 | 3.14163E+00 | 6.49921E-03 | +| 100000 | 3.14195E+00 | 4.11195E-01 | 3.14163E+00 | 6.36837E-03 | +| 1000000 | 3.14003E+00 | 4.14114E-01 | 3.14128E+00 | 6.39217E-03 | +| 10000000 | 3.14213E+00 | 4.13838E-01 | 3.14160E+00 | 6.40784E-03 | +|-----------------------------------------------------------------------| + +However, it is unfair to study one-dimensional integrals with MC methods! + +!eblock + + + + + +!split +===== Acceptance-Rejection Method ===== +!bblock +# add ref +This is a rather simple and appealing +method after von Neumann. Assume that we are looking at an interval +$x\in [a,b]$, this being the domain of the PDF $p(x)$. Suppose also that +the largest value our distribution function takes in this interval +is $M$, that is +!bt +\begin{equation*} + p(x) \le M \hspace{1cm} x\in [a,b]. +\end{equation*} +!et +Then we generate a random number $x$ from the uniform distribution +for $x\in [a,b]$ and a corresponding number $s$ for the uniform +distribution between $[0,M]$. If +!bt +\begin{equation*} +p(x) \ge s, +\end{equation*} +!et +we accept the new value of $x$, else we generate +again two new random numbers $x$ and $s$ and perform the test +in the latter equation again. +!eblock + +!split +===== Acceptance-Rejection Method ===== +!bblock + +As an example, consider the evaluation of the integral + +!bt +\begin{equation*} + I=\int_0^3\exp{(x)}dx. +\end{equation*} +!et +Obviously to derive a closed-form expression is much easier, however the integrand could pose some more +difficult challenges. The aim here is simply to show how to implent the acceptance-rejection algorithm. +The integral is the area below the curve $f(x)=\exp{(x)}$. If we uniformly fill the rectangle +spanned by $x\in [0,3]$ and $y\in [0,\exp{(3)}]$, the fraction below the curve obtained from a uniform distribution, and +multiplied by the area of the rectangle, should approximate the chosen integral. +!eblock + +!split +===== Acceptance-Rejection Method ===== +!bblock +It is rather +easy to implement this numerically, as shown in the following code. +!bc cppcod +// Loop over Monte Carlo trials n + integral =0.; + for ( int i = 1; i <= n; i++){ +// Finds a random value for x in the interval [0,3] + x = 3*ran0(&idum); +// Finds y-value between [0,exp(3)] + y = exp(3.0)*ran0(&idum); +// if the value of y at exp(x) is below the curve, we accept + if ( y < exp(x)) s = s+ 1.0; +// The integral is area enclosed below the line f(x)=exp(x) + } +// Then we multiply with the area of the rectangle and divide by the number of cycles + Integral = 3.*exp(3.)*s/n +!ec +!eblock + +!split +===== Monte Carlo Integration of Multidimensional Integrals ===== +!bblock +When we deal with multidimensional integrals of the form + +!bt +\begin{equation*} + I=\int_{a_1}^{b_1}dx_1\int_{a_2}^{b_2}dx_2\dots \int_{a_d}^{b_d}dx_d g(x_1,\dots,x_d), +\end{equation*} +!et +with +$x_i$ defined in the interval $[a_i,b_i]$ we would typically +need a transformation +of variables of the form + +!bt +\begin{equation*} + x_i=a_i+(b_i-a_i)t_i, +\end{equation*} +!et +if we were to use the uniform distribution on the interval $[0,1]$. +!eblock + +!split +===== Monte Carlo Integration of Multidimensional Integrals ===== +!bblock +In this case, we need a +Jacobi determinant + +!bt +\begin{equation*} + \prod_{i=1}^d (b_i-a_i), +\end{equation*} +!et +and to convert the function $g(x_1,\dots,x_d)$ to + +!bt +\begin{equation*} + g(x_1,\dots,x_d)\rightarrow + g(a_1+(b_1-a_1)t_1,\dots,a_d+(b_d-a_d)t_d). +\end{equation*} +!et +!eblock + +!split +===== Monte Carlo Integration of Multidimensional Integrals ===== +!bblock +As an example, consider the following six-dimensional +integral +!bt +\begin{equation*} + \int_{-\infty}^{\infty}\mathbf{dxdy}g(\mathbf{x, y}), +\end{equation*} +!et +where + +!bt +\begin{equation*} + g(\mathbf{x, y})=\exp{(-\mathbf{x}^2-\mathbf{y}^2)}(\mathbf{x}-\mathbf{y})^2 +\end{equation*} +!et +with $d=6$. +!eblock +!split +===== Monte Carlo Integration of Multidimensional Integrals ===== +!bblock +We can solve this integral by employing our brute force scheme, +or using importance sampling and random variables distributed +according to a gaussian PDF. For the latter, if we set +the mean value +$\mu=0$ and the standard deviation $\sigma=1/\sqrt{2}$, we have +!bt +\begin{equation*} + \frac{1}{\sqrt{\pi}}\exp{(-x^2)}, +\end{equation*} +!et +and using this normal distribution we rewrite our integral as +!bt +\begin{equation*} + \pi^3\int\prod_{i=1}^6\left( + \frac{1}{\sqrt{\pi}}\exp{(-x_i^2)}\right) + (\mathbf{x}-\mathbf{y})^2dx_1.\dots dx_6. +\end{equation*} +!et +!eblock + +!split +===== Monte Carlo Integration of Multidimensional Integrals ===== +!bblock + +We rewrite it in a more compact form as +!bt +\begin{equation*} + \int f(x_1,\dots,x_d)F(x_1,\dots,x_d)\prod_{i=1}^6dx_i, +\end{equation*} +!et +where $f$ is the above normal distribution +and +!bt +\begin{equation*} + F(x_1,\dots,x_6)=F(\mathbf{x, y})=(\mathbf{x}-\mathbf{y})^2, +\end{equation*} +!et +!eblock + +!split +===== Brute Force Integration ===== +!bblock +Below we list two codes, one for the brute force integration +and the other employing importance sampling with a gaussian distribution. +!bc cppcod +#include +#include +#include +#include "lib.h" +using namespace std; + +double brute_force_MC(double *); +// Main function begins here +int main() +{ + int n; + double x[6], y, fx; + double int_mc = 0.; double variance = 0.; + double sum_sigma= 0. ; long idum=-1 ; + double length = 5.; // we fix the max size of the box to L=5 + double jacobidet = pow((2*length),6); + cout << "Read in the number of Monte-Carlo samples" << endl; + cin >> n; +// evaluate the integral with importance sampling + for ( int i = 1; i <= n; i++){ +// x[] contains the random numbers for all dimensions + for (int j = 0; j< 6; j++) { + x[j]=-length+2*length*ran0(&idum); + } + fx=brute_force_MC(x); + int_mc += fx; + sum_sigma += fx*fx; + } + int_mc = int_mc/((double) n ); + sum_sigma = sum_sigma/((double) n ); + variance=sum_sigma-int_mc*int_mc; +// final output + cout << setiosflags(ios::showpoint | ios::uppercase); + cout << " Monte carlo result= " << setw(10) << setprecision(8) << jacobidet*int_mc; + cout << " Sigma= " << setw(10) << setprecision(8) << volume*sqrt(variance/((double) n )) << endl; + return 0; +} // end of main program + +// this function defines the integrand to integrate + +double brute_force_MC(double *x) +{ +// evaluate the different terms of the exponential + double xx=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]; + double yy=x[3]*x[3]+x[4]*x[4]+x[5]*x[5]; + double xy=pow((x[0]-x[3]),2)+pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2); + return exp(-xx-yy)*xy; +} // end function for the integrand +!ec +!eblock + +!split +===== Importance Sampling ===== +!bblock +This code includes a call to the function $normal_random$, which produces +random numbers from a gaussian distribution. +!bc cppcod +// importance sampling with gaussian deviates +#include +#include +#include +#include "lib.h" +using namespace std; + +double gaussian_MC(double *); +double gaussian_deviate(long *); +// Main function begins here +int main() +{ + int n; + double x[6], y, fx; + cout << "Read in the number of Monte-Carlo samples" << endl; + cin >> n; + double int_mc = 0.; double variance = 0.; + double sum_sigma= 0. ; long idum=-1 ; + double jacobidet = pow(acos(-1.),3.); + double sqrt2 = 1./sqrt(2.); +// evaluate the integral with importance sampling + for ( int i = 1; i <= n; i++){ +// x[] contains the random numbers for all dimensions + for (int j = 0; j < 6; j++) { + x[j] = gaussian_deviate(&idum)*sqrt2; + } + fx=gaussian_MC(x); + int_mc += fx; + sum_sigma += fx*fx; + } + int_mc = int_mc/((double) n ); + sum_sigma = sum_sigma/((double) n ); + variance=sum_sigma-int_mc*int_mc; +// final output + cout << setiosflags(ios::showpoint | ios::uppercase); + cout << " Monte carlo result= " << setw(10) << setprecision(8) << jacobidet*int_mc; + cout << " Sigma= " << setw(10) << setprecision(8) << volume*sqrt(variance/((double) n )) << endl; + return 0; +} // end of main program + +// this function defines the integrand to integrate + +double gaussian_MC(double *x) +{ +// evaluate the different terms of the exponential + double xy=pow((x[0]-x[3]),2)+pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2); + return xy; +} // end function for the integrand + +// random numbers with gaussian distribution +double gaussian_deviate(long * idum) +{ + static int iset = 0; + static double gset; + double fac, rsq, v1, v2; + + if ( idum < 0) iset =0; + if (iset == 0) { + do { + v1 = 2.*ran0(idum) -1.0; + v2 = 2.*ran0(idum) -1.0; + rsq = v1*v1+v2*v2; + } while (rsq >= 1.0 || rsq == 0.); + fac = sqrt(-2.*log(rsq)/rsq); + gset = v1*fac; + iset = 1; + return v2*fac; + } else { + iset =0; + return gset; + } +} // end function for gaussian deviates +!ec +!eblock + + +!split +===== Python codes ===== +!bblock +The first code here is an example of a python which computes the above integral using +the brute force approach +!bc pyscpro +#Monte Carlo integration in 6 dimensions + +import numpy,math +import sys + +def integrand(x): + """Calculates the integrand + exp(-a*(x1^2+x2^2+...+x6^2)-b*[(x1-x4)^2+...+(x3-x6)^2]) + from the values in the 6-dimensional array x.""" + a = 1.0 + b = 0.5 + + x2 = numpy.sum(x**2) + xy = (x[0]-x[3])**2 + (x[1]-x[4])**2 + (x[2]-x[5])**2 + + return numpy.exp(-a*x2-b*xy) + +#Main program + +#Integration limits: x[i] in (-5, 5) +L = 5.0 +jacobi = (2*L)**6 + +N = 100000 + +#Evaluate the integral +sum = 0.0 +sum2 = 0.0 +for i in xrange(N): + #Generate random coordinates to sample at + x = numpy.array([-L+2*L*numpy.random.random() for j in xrange(6)]) + + fx = integrand(x) + sum += fx + sum2 += fx**2 +#Calculate expt. values for fx and fx^2 +sum /=float(N) +sum2/=float(N) + +#Result +int_mc = jacobi*sum; +#Gaussian standard deviation +sigma = jacobi*math.sqrt((sum2-sum**2)/float(N)) + +#Output +print "Montecarlo result = %10.8E" % int_mc +print "Sigma = %10.8E" % sigma + +!ec +!eblock + +!split +===== Python codes, importance sampling ===== +!bblock +The second code, displayed here, uses importance sampling and random numbers that follow the normal distribution +the brute force approach +!bc pyscpro +#Monte Carlo integration with importance sampling + +import numpy,math +import sys + +def integrand(x): + """Calculates the integrand + exp(-b*[(x1-x4)^2+...+(x3-x6)^2]) + from the values in the 6-dimensional array x.""" + b = 0.5 + + xy = (x[0]-x[3])**2 + (x[1]-x[4])**2 + (x[2]-x[5])**2 + return numpy.exp(-b*xy) + +#Main program + +#Jacobi determinant +jacobi = math.acos(-1.0)**3 +sqrt2 = 1.0/math.sqrt(2) + +N = 100000 + +#Evaluate the integral +sum = 0.0 +sum2 = 0.0 +for i in xrange(N): + #Generate random, gaussian distributed coordinates to sample at + x = numpy.array([numpy.random.normal()*sqrt2 for j in xrange(6)]) + + fx = integrand(x) + sum += fx + sum2 += fx**2 +#Calculate expt. values for fx and fx^2 +sum /=float(N) +sum2/=float(N) + +#Result +int_mc = jacobi*sum; +#Gaussian standard deviation +sigma = jacobi*math.sqrt((sum2-sum**2)/float(N)) + +#Output +print "Montecarlo result = %10.8E" % int_mc +print "Sigma = %10.8E" % sigma +!ec +!eblock + + + diff --git a/doc/src/Statistics/clean.sh b/doc/src/Statistics/clean.sh new file mode 100644 index 000000000..2e5da2c72 --- /dev/null +++ b/doc/src/Statistics/clean.sh @@ -0,0 +1,3 @@ +#!/bin/sh +doconce clean +rm -rf *.pdf *.tex ipynb*.tar.gz *.html ._*.html *~ reveal.js Trash README.txt diff --git a/doc/src/Statistics/datafiles/autocor.dat b/doc/src/Statistics/datafiles/autocor.dat new file mode 100644 index 000000000..8ee8c5fe3 --- /dev/null +++ b/doc/src/Statistics/datafiles/autocor.dat @@ -0,0 +1,1000 @@ + 0 1.0000000 + 1 0.030436777 + 2 -0.041967823 + 3 0.0025843740 + 4 -0.013971981 + 5 0.031863200 + 6 0.011244624 + 7 -0.0017268302 + 8 -0.032367877 + 9 0.0048970777 + 10 -0.011108971 + 11 0.029249261 + 12 -0.060102972 + 13 -0.010339285 + 14 -0.016163202 + 15 -0.014328605 + 16 0.037522532 + 17 0.0086741695 + 18 0.042486236 + 19 0.038599298 + 20 -0.017576480 + 21 0.026784371 + 22 -0.068080722 + 23 -0.032413624 + 24 -0.033787280 + 25 0.0044378763 + 26 0.0084675258 + 27 -0.049112322 + 28 0.029052476 + 29 -0.010253671 + 30 -0.020415539 + 31 -0.026030116 + 32 0.035646939 + 33 -0.016364511 + 34 0.031223752 + 35 -0.00098907934 + 36 0.011793763 + 37 0.044375824 + 38 -0.0075893899 + 39 0.038505140 + 40 -0.050820678 + 41 -0.0058587490 + 42 -0.014352719 + 43 -0.030117666 + 44 -0.010670435 + 45 0.037061620 + 46 0.012236259 + 47 0.012275474 + 48 -0.043753664 + 49 -0.046814374 + 50 -0.0078373590 + 51 0.030224467 + 52 -0.0014576041 + 53 -0.0073935525 + 54 -0.0015361878 + 55 0.0063681276 + 56 0.00077378384 + 57 0.027396611 + 58 0.026288042 + 59 0.041165848 + 60 0.033625798 + 61 -0.021151787 + 62 0.0091334499 + 63 0.019146160 + 64 -0.010915406 + 65 -0.038860498 + 66 -0.0086368193 + 67 0.0021648597 + 68 0.026319537 + 69 0.0053794411 + 70 0.039529055 + 71 0.0087119709 + 72 -0.031583705 + 73 0.040882208 + 74 0.0081181815 + 75 0.0016127077 + 76 0.057618438 + 77 -0.020863018 + 78 -0.053206466 + 79 -0.048595441 + 80 -0.0019519214 + 81 -0.0091520929 + 82 -0.0093091454 + 83 -0.025526629 + 84 0.0038194868 + 85 0.0056035808 + 86 -0.041792435 + 87 0.019370857 + 88 0.014560618 + 89 -0.016079104 + 90 0.017263262 + 91 0.030090005 + 92 0.046198390 + 93 0.014260591 + 94 0.0035178537 + 95 0.014851525 + 96 -0.015769457 + 97 -0.0019402466 + 98 0.022252302 + 99 -0.017622777 + 100 -0.029172137 + 101 0.00092111483 + 102 -0.019516339 + 103 -0.014110995 + 104 -0.0053766072 + 105 0.015551435 + 106 -0.0092984599 + 107 -0.017125844 + 108 -0.023354720 + 109 0.037924881 + 110 -0.013205795 + 111 -0.073863600 + 112 -0.014986822 + 113 0.022838581 + 114 0.028783448 + 115 -0.0046507703 + 116 0.0072623794 + 117 -0.017007609 + 118 0.0054601073 + 119 0.014645594 + 120 -0.018330664 + 121 -0.0074984609 + 122 -0.019475491 + 123 0.041232837 + 124 0.035696824 + 125 -0.010334698 + 126 -0.011736466 + 127 0.030226422 + 128 0.021398621 + 129 0.024604001 + 130 -0.016037954 + 131 0.038544269 + 132 0.061300628 + 133 0.019470825 + 134 -0.011061461 + 135 0.049068767 + 136 -0.049298447 + 137 0.0039132268 + 138 0.020882292 + 139 -0.025826097 + 140 -0.011796518 + 141 0.028217876 + 142 0.016595503 + 143 0.025935569 + 144 -0.0078950322 + 145 0.0010445150 + 146 -0.019355185 + 147 -0.026751109 + 148 0.023582576 + 149 0.025495541 + 150 -0.0046856525 + 151 -0.0053126503 + 152 -0.018928898 + 153 0.017209905 + 154 -0.020543359 + 155 -0.0019265343 + 156 -0.0045039479 + 157 -0.044022012 + 158 -0.0054810482 + 159 0.0097895871 + 160 0.069425506 + 161 0.0074504423 + 162 -0.016934482 + 163 -0.0075259798 + 164 -0.00021461104 + 165 -0.056810813 + 166 0.0041571957 + 167 -0.051444304 + 168 0.018242092 + 169 -0.0092588016 + 170 -0.058149110 + 171 -0.019497092 + 172 0.033414724 + 173 0.019228122 + 174 -0.017014927 + 175 0.023859853 + 176 0.0093103834 + 177 0.033798290 + 178 0.013270749 + 179 0.067569382 + 180 -0.010644310 + 181 -0.046364881 + 182 0.018396640 + 183 0.020785098 + 184 -0.021406123 + 185 -0.022586765 + 186 -0.015545857 + 187 -0.030076727 + 188 -0.0065117780 + 189 -0.0040088089 + 190 0.010598275 + 191 -0.017176792 + 192 0.056184566 + 193 -0.011321257 + 194 -0.0029512385 + 195 -0.014918170 + 196 -0.0039190314 + 197 0.056045986 + 198 -0.032240800 + 199 -0.024267093 + 200 0.025936442 + 201 0.024728180 + 202 -0.00089880027 + 203 -0.029831804 + 204 -0.027023166 + 205 -0.033939923 + 206 -0.033440421 + 207 0.026431974 + 208 0.021336927 + 209 -0.052701327 + 210 -0.036889376 + 211 0.050380536 + 212 -0.011311638 + 213 0.0057673144 + 214 -0.0027656116 + 215 0.042261502 + 216 0.028936072 + 217 0.028579634 + 218 0.024570794 + 219 -0.0069215979 + 220 -0.019236105 + 221 -0.021460579 + 222 0.045327931 + 223 0.017796223 + 224 -0.021253081 + 225 0.0071524123 + 226 -0.015126679 + 227 -0.051499638 + 228 0.012432292 + 229 0.0074777375 + 230 -0.017758289 + 231 0.014844143 + 232 0.013394570 + 233 -0.018637091 + 234 -0.0099195677 + 235 -0.011840993 + 236 0.056187841 + 237 0.015187364 + 238 -0.027196415 + 239 -0.021471154 + 240 -0.015489221 + 241 0.0011648445 + 242 -0.014737635 + 243 0.041224909 + 244 -0.048386598 + 245 -0.028771515 + 246 -0.0069030430 + 247 -0.0082757555 + 248 0.035867587 + 249 -0.028705115 + 250 -0.023521399 + 251 -0.0080840144 + 252 0.032038834 + 253 0.0016768482 + 254 0.017568742 + 255 -0.038854623 + 256 -0.019572338 + 257 -0.018945188 + 258 -0.044047531 + 259 0.035712621 + 260 0.0052305122 + 261 -0.019221249 + 262 0.016556546 + 263 -0.018528132 + 264 -0.033294901 + 265 -0.029063893 + 266 0.068630731 + 267 0.014630449 + 268 -0.010624595 + 269 0.0068084665 + 270 -0.013154851 + 271 -0.0031250199 + 272 -0.022697076 + 273 0.052943127 + 274 0.013646959 + 275 0.019661200 + 276 -0.0021662857 + 277 0.0048494673 + 278 0.0093039642 + 279 -0.050613478 + 280 0.0037930281 + 281 0.0027056109 + 282 0.0067641602 + 283 -0.031584591 + 284 -0.016064322 + 285 -0.0010516554 + 286 -0.028260447 + 287 0.0040810671 + 288 0.010358807 + 289 -0.019334825 + 290 -0.048260958 + 291 0.035049342 + 292 0.040489272 + 293 0.0020011032 + 294 -0.047585009 + 295 -0.0069058676 + 296 -0.0068136076 + 297 -0.025701037 + 298 0.016612923 + 299 0.00069465478 + 300 -0.025587862 + 301 0.029814669 + 302 -0.041293387 + 303 -0.017292295 + 304 -0.013939533 + 305 -0.017807378 + 306 0.035623264 + 307 0.023870958 + 308 -0.0048061196 + 309 0.013506371 + 310 0.024774288 + 311 0.0073845590 + 312 0.0074994639 + 313 0.011664191 + 314 0.039796920 + 315 -0.021968341 + 316 0.030270961 + 317 0.011032615 + 318 -0.013020314 + 319 -0.016347247 + 320 -0.0053552657 + 321 -0.0018146771 + 322 -0.028324460 + 323 -0.033892233 + 324 0.040970876 + 325 0.013724516 + 326 -0.025758979 + 327 0.019040388 + 328 -0.0047334635 + 329 -0.0034112772 + 330 0.0068517387 + 331 -0.026995154 + 332 -0.014256522 + 333 0.019176374 + 334 0.0051692845 + 335 0.041770495 + 336 -0.0019396636 + 337 0.010891816 + 338 -0.044531596 + 339 -0.0077558642 + 340 -0.025102912 + 341 -0.010674026 + 342 0.024980787 + 343 0.0024171666 + 344 -0.0081651699 + 345 0.023649184 + 346 0.0069294682 + 347 -0.039211268 + 348 0.021395177 + 349 0.010078175 + 350 0.0069000193 + 351 0.028906152 + 352 0.022775633 + 353 -0.023275057 + 354 0.017183626 + 355 -0.0034314696 + 356 0.010874216 + 357 -0.029315379 + 358 0.010861994 + 359 0.017521148 + 360 0.025805095 + 361 0.019252809 + 362 -0.064531812 + 363 0.048847457 + 364 0.033641876 + 365 -0.0051765311 + 366 0.0091743379 + 367 0.033470461 + 368 0.016391479 + 369 0.033188930 + 370 0.00020395242 + 371 0.010477670 + 372 0.043657980 + 373 -0.045721493 + 374 0.0095296094 + 375 0.038701345 + 376 -0.019764861 + 377 -0.0024760707 + 378 -0.035616363 + 379 -0.025423739 + 380 0.0048029114 + 381 -0.0057620708 + 382 -0.015241488 + 383 -0.0094964164 + 384 -0.015301911 + 385 -0.029335964 + 386 -0.0085647523 + 387 -0.035018913 + 388 -0.018494103 + 389 -0.025290309 + 390 0.043807258 + 391 0.024917719 + 392 -0.037659446 + 393 -0.0094455059 + 394 -0.0092903549 + 395 0.021372140 + 396 0.0029776435 + 397 0.048249941 + 398 -0.051252331 + 399 -0.042462183 + 400 -0.027531735 + 401 0.032780000 + 402 -0.012075698 + 403 0.017104870 + 404 0.0056918634 + 405 -0.037383540 + 406 -0.0099994917 + 407 0.034147318 + 408 0.035420675 + 409 -0.018830626 + 410 0.0057399934 + 411 0.028717444 + 412 0.0043248247 + 413 -0.015179249 + 414 0.038564990 + 415 0.028899559 + 416 0.0070276321 + 417 -0.020818426 + 418 -0.0013244685 + 419 -0.016868518 + 420 0.015794320 + 421 -0.034150397 + 422 -0.011318810 + 423 0.0093680067 + 424 -0.0056974750 + 425 0.00089406498 + 426 0.023839990 + 427 -0.023249489 + 428 0.0036581356 + 429 -0.0068173069 + 430 -0.010542605 + 431 0.021106295 + 432 -0.026663001 + 433 -0.00065540509 + 434 0.0056271702 + 435 0.014968416 + 436 0.0069719584 + 437 -0.032506238 + 438 -0.018610431 + 439 0.013945493 + 440 0.030296160 + 441 0.015990399 + 442 0.027637874 + 443 -0.040348487 + 444 -0.0082368848 + 445 0.0084959142 + 446 -0.023774043 + 447 -0.021734126 + 448 0.012168605 + 449 -0.0074461207 + 450 -0.016125466 + 451 0.033331887 + 452 0.00062736237 + 453 -0.012717250 + 454 0.017806603 + 455 -0.027145155 + 456 0.0012005238 + 457 -0.0077436248 + 458 -0.028097243 + 459 -0.00095474376 + 460 -0.0014605876 + 461 -0.035595271 + 462 -0.051237530 + 463 0.010146341 + 464 -0.028954767 + 465 0.0037797756 + 466 -0.030882319 + 467 0.047912444 + 468 0.011363212 + 469 -0.0085205849 + 470 -0.0039355171 + 471 -0.031696814 + 472 0.012773460 + 473 0.0055999575 + 474 -0.00037892267 + 475 0.015504229 + 476 -0.0076580626 + 477 0.045706707 + 478 0.0071347194 + 479 -0.045884838 + 480 -0.043718227 + 481 0.018240807 + 482 0.054386867 + 483 0.046911775 + 484 0.045204253 + 485 0.010090395 + 486 0.026767753 + 487 0.010723149 + 488 0.025154764 + 489 -0.024373964 + 490 -0.018122001 + 491 0.018975902 + 492 -0.029915565 + 493 0.021781915 + 494 -0.012285808 + 495 0.0027546898 + 496 0.00062612064 + 497 -0.040850871 + 498 -0.023513739 + 499 0.021760239 + 500 -0.021001688 + 501 0.0096979566 + 502 0.0063539835 + 503 0.021533295 + 504 0.023511110 + 505 -0.012760504 + 506 -0.018141212 + 507 -0.021346409 + 508 -0.0022698665 + 509 0.019874015 + 510 0.011557117 + 511 -0.042294961 + 512 0.029066005 + 513 0.0017714648 + 514 3.3189262E-05 + 515 -0.013922936 + 516 -0.0045819016 + 517 -0.023030698 + 518 0.017014357 + 519 -0.015759868 + 520 0.022547087 + 521 0.010541576 + 522 -0.023908324 + 523 -0.012544793 + 524 -0.034918466 + 525 -0.032696570 + 526 -0.014898371 + 527 0.024734697 + 528 0.029617332 + 529 0.0092461414 + 530 -0.035946729 + 531 -0.011867012 + 532 0.034753028 + 533 3.7608681E-05 + 534 -0.026655371 + 535 -0.0085675831 + 536 -0.0043169745 + 537 0.00019284442 + 538 -0.0027829108 + 539 0.014288886 + 540 0.0036864577 + 541 -0.016568991 + 542 0.011661583 + 543 0.010312897 + 544 0.047191126 + 545 0.010246982 + 546 0.018691967 + 547 0.012782215 + 548 0.010575964 + 549 0.043346000 + 550 0.0012075999 + 551 0.0048928567 + 552 0.0088800511 + 553 -0.026749858 + 554 -0.049468570 + 555 0.0087226114 + 556 -0.025355412 + 557 0.00041911744 + 558 -0.0039277925 + 559 -0.019867410 + 560 0.0090696073 + 561 -0.020407695 + 562 -0.028022150 + 563 0.00021588069 + 564 0.0036235564 + 565 0.012238820 + 566 0.0016474398 + 567 0.019772585 + 568 -0.0079130745 + 569 -0.0048698647 + 570 -0.00046839106 + 571 -0.0068985684 + 572 0.0038773945 + 573 0.025536952 + 574 0.014407956 + 575 -0.016018662 + 576 0.018400000 + 577 -0.0074673768 + 578 0.014582128 + 579 0.015277437 + 580 -0.010662929 + 581 -0.028007675 + 582 -0.025491028 + 583 0.010355834 + 584 -0.017664380 + 585 -0.0027916888 + 586 -0.0031681423 + 587 -0.020491989 + 588 0.030287661 + 589 -0.0018555509 + 590 -0.00029702455 + 591 0.0077901131 + 592 -0.0074424409 + 593 -0.042697622 + 594 -0.0081579992 + 595 0.0075807326 + 596 0.0062213137 + 597 -0.024978210 + 598 -0.015466545 + 599 -0.021531833 + 600 -0.014845678 + 601 -0.0047515868 + 602 0.026920902 + 603 0.0093622167 + 604 0.035153363 + 605 0.0062605573 + 606 -0.0025669440 + 607 -0.019158868 + 608 0.033285966 + 609 -0.0054035913 + 610 -0.017200913 + 611 0.014100356 + 612 -0.0090120505 + 613 0.0039391323 + 614 -0.028582199 + 615 0.026504777 + 616 -0.031792297 + 617 0.0053277806 + 618 0.010772900 + 619 -0.013054304 + 620 0.012242009 + 621 -0.065873242 + 622 0.0032291061 + 623 0.0091991992 + 624 0.00032081972 + 625 0.033746878 + 626 0.021605607 + 627 -0.00065569553 + 628 -0.0052540101 + 629 -0.0088500047 + 630 0.016983225 + 631 0.017762578 + 632 -0.020731913 + 633 0.027829821 + 634 0.011591300 + 635 -0.028061398 + 636 0.034537856 + 637 -0.0053690020 + 638 -0.045885706 + 639 -0.017892409 + 640 -0.010396254 + 641 0.028504922 + 642 0.024906302 + 643 0.031601151 + 644 -0.020471528 + 645 0.012436367 + 646 -0.0060652007 + 647 0.016340730 + 648 -0.051146749 + 649 0.025541636 + 650 -0.0091187569 + 651 -0.028358025 + 652 0.0090547667 + 653 -0.039390098 + 654 -0.025279051 + 655 -0.0098439764 + 656 0.010135984 + 657 -0.0099451758 + 658 -0.013868539 + 659 0.011209477 + 660 -0.015549717 + 661 0.015580477 + 662 0.034790993 + 663 0.036386181 + 664 0.017019255 + 665 -0.0071014117 + 666 -0.026363213 + 667 -0.0011162350 + 668 0.010430304 + 669 -0.0037094956 + 670 -0.017952290 + 671 -0.0077706804 + 672 0.0029775795 + 673 -0.0061743735 + 674 -0.015180197 + 675 -0.0090826895 + 676 0.013756556 + 677 0.025269597 + 678 0.034896932 + 679 0.0045113518 + 680 -0.0027466582 + 681 -0.0020897520 + 682 0.0077071015 + 683 0.017099343 + 684 -0.0015647601 + 685 -0.042943759 + 686 -0.019877423 + 687 -0.014563904 + 688 0.028812347 + 689 0.013691120 + 690 -0.024519950 + 691 -0.0044018710 + 692 0.0078443611 + 693 0.027075429 + 694 -0.012549527 + 695 0.016024749 + 696 -0.0046314243 + 697 0.0084472590 + 698 -0.015057447 + 699 0.026575886 + 700 0.010553947 + 701 -0.0033593037 + 702 0.0030676930 + 703 0.0086218092 + 704 0.0037965560 + 705 -0.014944414 + 706 -0.0029175505 + 707 -0.022225334 + 708 0.0035673060 + 709 -0.014910503 + 710 -0.0028479147 + 711 -0.0012980229 + 712 -0.024413269 + 713 -0.014945607 + 714 -0.00092249340 + 715 0.0060404214 + 716 -0.018852282 + 717 -0.0093385434 + 718 0.0040141908 + 719 0.0024247280 + 720 0.018611644 + 721 -0.0052456314 + 722 -0.022206528 + 723 0.012921902 + 724 0.00064998606 + 725 0.0043656363 + 726 -0.0026142811 + 727 0.0066569906 + 728 -0.014715930 + 729 0.023061889 + 730 -0.014595305 + 731 0.018909306 + 732 0.0019405656 + 733 -0.024034112 + 734 0.044901386 + 735 -0.012970072 + 736 0.011127732 + 737 -0.014563174 + 738 -0.0064960244 + 739 0.0023165471 + 740 -0.00087576396 + 741 -0.0049030238 + 742 -0.019033070 + 743 -0.0021911524 + 744 0.025625643 + 745 -0.014500976 + 746 -0.012423680 + 747 -0.017058080 + 748 0.0045641993 + 749 0.0084538383 + 750 0.0021116607 + 751 -0.0033408150 + 752 -0.0072360325 + 753 -0.0022798159 + 754 -0.021979932 + 755 0.023315014 + 756 0.0052468019 + 757 0.016949879 + 758 -0.017138783 + 759 0.0078100113 + 760 0.0099470522 + 761 0.0096305195 + 762 0.027398820 + 763 -0.025646831 + 764 0.013123399 + 765 0.015549034 + 766 0.0023558555 + 767 -0.0096405092 + 768 0.012253011 + 769 -0.0051483047 + 770 -0.0084922180 + 771 -0.012338577 + 772 -0.012115014 + 773 -0.0054183951 + 774 0.015587209 + 775 0.015120992 + 776 -0.012270671 + 777 -0.0034455014 + 778 0.010222090 + 779 -0.020718897 + 780 -0.033416426 + 781 0.031396462 + 782 -0.016828081 + 783 0.0037494749 + 784 -0.0072224179 + 785 0.0056141543 + 786 0.018963466 + 787 -0.0082246575 + 788 -0.0050421303 + 789 -0.022400808 + 790 0.011845444 + 791 0.0037973135 + 792 0.0081527811 + 793 0.0013800440 + 794 0.019983158 + 795 -0.0056201910 + 796 -0.0043737064 + 797 -0.0091988539 + 798 -0.0098326674 + 799 -0.024412160 + 800 -0.0057559680 + 801 0.0066857351 + 802 0.021505241 + 803 0.0030792036 + 804 -0.019226409 + 805 0.0045538061 + 806 -0.0056070254 + 807 -0.0016601038 + 808 -0.024716882 + 809 0.0054226135 + 810 -0.00025770775 + 811 0.015113434 + 812 -0.0032290640 + 813 -0.027790661 + 814 -0.0041114266 + 815 0.013470213 + 816 -0.0092058822 + 817 -0.0095842170 + 818 0.014469428 + 819 0.019165159 + 820 0.0028352364 + 821 -9.5780268E-05 + 822 -0.0076164534 + 823 -0.013116541 + 824 0.011262012 + 825 -0.014590759 + 826 0.0072059449 + 827 -0.015742103 + 828 0.013380709 + 829 -0.027024814 + 830 -0.011718400 + 831 0.0070073072 + 832 -0.0058391392 + 833 -0.0041303314 + 834 -0.031348383 + 835 -0.010022600 + 836 -0.0021030126 + 837 -0.0088289877 + 838 0.011174023 + 839 0.0078621381 + 840 -0.020319884 + 841 -0.0054612370 + 842 0.0036168037 + 843 -0.021466662 + 844 0.017063947 + 845 -0.015876318 + 846 -0.0015561464 + 847 0.018739482 + 848 0.0070220237 + 849 0.010466066 + 850 -0.0065905279 + 851 0.0050221836 + 852 -0.0049380275 + 853 -0.0022540120 + 854 -0.0066911228 + 855 -0.00085489296 + 856 0.0091311886 + 857 0.0087961910 + 858 0.0084225981 + 859 0.0033897841 + 860 -0.0084094261 + 861 -0.0066896747 + 862 0.0077203999 + 863 0.0012605110 + 864 -0.0026422053 + 865 0.0036337117 + 866 0.0017212519 + 867 -0.015745443 + 868 -0.0092171165 + 869 -0.00011985970 + 870 -0.021726538 + 871 -0.0073536601 + 872 7.2163303E-05 + 873 0.0011758806 + 874 0.014912380 + 875 0.0012475430 + 876 0.022576078 + 877 0.0045818387 + 878 0.010462266 + 879 0.025107039 + 880 -0.013052848 + 881 0.0025105745 + 882 -0.010884996 + 883 -0.0025679305 + 884 -0.00040656738 + 885 0.014211879 + 886 -0.013404473 + 887 -0.014505338 + 888 -0.0027969620 + 889 0.0092900832 + 890 0.012893798 + 891 -0.011269371 + 892 0.016874925 + 893 -0.0047647532 + 894 0.0081446665 + 895 0.0058947644 + 896 0.00076487417 + 897 -0.0025185792 + 898 -0.013129512 + 899 -0.0067449885 + 900 -0.014411868 + 901 -0.011283279 + 902 0.011162838 + 903 -0.010418718 + 904 0.0032281419 + 905 -0.010894860 + 906 0.0011494465 + 907 0.0084593562 + 908 0.0065085104 + 909 0.014711201 + 910 0.0054832203 + 911 0.0012882653 + 912 -0.00089841363 + 913 0.0050714575 + 914 -0.011098348 + 915 0.0055713307 + 916 0.0024382259 + 917 -0.0074570407 + 918 -0.0020684133 + 919 0.00048080589 + 920 0.0021420954 + 921 0.0054954627 + 922 0.0021075997 + 923 0.0033861839 + 924 -0.015575036 + 925 -0.00010693208 + 926 0.0076125254 + 927 0.0033302873 + 928 -0.013058275 + 929 0.0013763743 + 930 -0.014574421 + 931 -0.0043236772 + 932 0.0015626567 + 933 -0.0033933812 + 934 0.010242569 + 935 -0.0067711917 + 936 -0.00072185958 + 937 0.0081792143 + 938 -0.0024182634 + 939 0.0038473932 + 940 0.0069120137 + 941 -0.00085246271 + 942 -0.0067411485 + 943 0.0029968194 + 944 0.0039549069 + 945 -0.0037928923 + 946 -0.0036476329 + 947 0.0040189508 + 948 -0.012532289 + 949 0.0058993881 + 950 -0.0053575605 + 951 0.0060915122 + 952 0.0013451894 + 953 -0.0045829537 + 954 0.018097604 + 955 -0.012682236 + 956 0.0020792812 + 957 0.0044188576 + 958 -0.012281183 + 959 0.0011869564 + 960 0.0043558809 + 961 -0.0042990975 + 962 -0.0068498043 + 963 0.014101593 + 964 -0.011923454 + 965 -0.00064539439 + 966 -0.0039817352 + 967 0.00014351008 + 968 0.0051525695 + 969 -0.00033345980 + 970 0.015672329 + 971 -0.0031266958 + 972 0.00093231204 + 973 -0.0034217829 + 974 -0.0027825064 + 975 -0.0058773123 + 976 0.0039172115 + 977 0.0022204348 + 978 -0.0019712632 + 979 0.0074446894 + 980 -0.0050907804 + 981 0.0058000460 + 982 -0.0064782460 + 983 -0.00086478009 + 984 0.0028678024 + 985 -0.0047832128 + 986 0.0071608110 + 987 0.0022478887 + 988 -0.0024566931 + 989 0.00024711205 + 990 -0.00029949750 + 991 0.00033041933 + 992 -0.0015456286 + 993 -0.0023284185 + 994 -0.0012638969 + 995 -0.00034855887 + 996 -0.0013436311 + 997 0.0024342189 + 998 0.0013168498 + 999 0.0019766962 diff --git a/doc/src/Statistics/datafiles/autocorrelation.cpp b/doc/src/Statistics/datafiles/autocorrelation.cpp new file mode 100644 index 000000000..bd11734ac --- /dev/null +++ b/doc/src/Statistics/datafiles/autocorrelation.cpp @@ -0,0 +1,67 @@ +// This function computes the autocorrelation function for +// the standard c++ random number generator + +#include +#include +#include +#include +using namespace std; +// output file as global variable +ofstream ofile; + +// Main function begins here +int main(int argc, char* argv[]) +{ + int n; + char *outfilename; + + cin >> n; + double MCint = 0.; double MCintsqr2=0.; + double invers_period = 1./RAND_MAX; // initialise the random number generator + srand(time(NULL)); // This produces the so-called seed in MC jargon + // Compute the variance and the mean value of the uniform distribution + // Compute also the specific values x for each cycle in order to be able to + // the covariance and the correlation function + // Read in output file, abort if there are too few command-line arguments + if( argc <= 2 ){ + cout << "Bad Usage: " << argv[0] << + " read also output file and number of cycles on same line" << endl; + exit(1); + } + else{ + outfilename=argv[1]; + } + ofile.open(outfilename); + // Get the number of Monte-Carlo samples + n = atoi(argv[2]); + double *X; + X = new double[n]; + for (int i = 0; i < n; i++){ + double x = double(rand())*invers_period; + X[i] = x; + MCint += x; + MCintsqr2 += x*x; + } + double Mean = MCint/((double) n ); + MCintsqr2 = MCintsqr2/((double) n ); + double STDev = sqrt(MCintsqr2-Mean*Mean); + double Variance = MCintsqr2-Mean*Mean; +// Write mean value and standard deviation + cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl; + + // Now we compute the autocorrelation function, setting the distance d between two + // to a most 1/4 of the total number of cycles + double *autocor; autocor = new double[n]; + for (int j = 0; j < n; j++){ + double sum = 0.0; + for (int k = 0; k < (n-j); k++){ + sum += (X[k]-Mean)*(X[k+j]-Mean); + } + autocor[j] = sum/Variance/((double) n ); + ofile << setiosflags(ios::showpoint | ios::uppercase); + ofile << setw(15) << setprecision(8) << j; + ofile << setw(15) << setprecision(8) << autocor[j] << endl; + } + ofile.close(); // close output file + return 0; +} // end of main program diff --git a/doc/src/Statistics/datafiles/automersenne.cpp b/doc/src/Statistics/datafiles/automersenne.cpp new file mode 100644 index 000000000..f4039a3b9 --- /dev/null +++ b/doc/src/Statistics/datafiles/automersenne.cpp @@ -0,0 +1,76 @@ +// This function computes the autocorrelation function for +// the standard c++ random number generator + +#include +#include +#include +#include +#include + +using namespace std; +// output file as global variable +ofstream ofile; + +// Main function begins here +int main(int argc, char* argv[]) +{ + int n; + char *outfilename; + + cin >> n; + double MCint = 0.; double MCintsqr2=0.; + // Initialize the seed and call the Mersienne algo + std::random_device rd; + std::mt19937_64 gen(rd()); + // Set up the uniform distribution for x \in [[0, 1] + std::uniform_real_distribution RandomNumberGenerator(0.0,1.0); + // Compute the variance and the mean value of the uniform distribution + // Compute also the specific values x for each cycle in order to be able to + // the covariance and the correlation function + // Read in output file, abort if there are too few command-line arguments + if( argc <= 2 ){ + cout << "Bad Usage: " << argv[0] << + " read also output file and number of cycles on same line" << endl; + exit(1); + } + else{ + outfilename=argv[1]; + } + ofile.open(outfilename); + // Get the number of Monte-Carlo samples + n = atoi(argv[2]); + double *X; + X = new double[n]; + for (int i = 0; i < n; i++){ + double x = RandomNumberGenerator(gen); + X[i] = x; + MCint += x; + MCintsqr2 += x*x; + } + double Mean = MCint/((double) n ); + MCintsqr2 = MCintsqr2/((double) n ); + double STDev = sqrt(MCintsqr2-Mean*Mean); + double Variance = MCintsqr2-Mean*Mean; +// Write mean value and standard deviation + cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl; + + // Now we compute the autocorrelation function, setting the distance d + double *autocor; autocor = new double[n]; + for (int j = 0; j < n; j++){ + double sum = 0.0; + for (int k = 0; k < (n-j); k++){ + sum += (X[k]-Mean)*(X[k+j]-Mean); + } + autocor[j] = sum/Variance/((double) n ); + ofile << setiosflags(ios::showpoint | ios::uppercase); + ofile << setw(15) << setprecision(8) << j; + ofile << setw(15) << setprecision(8) << autocor[j] << endl; + } + ofile.close(); // close output file + return 0; +} // end of main program + + + + + diff --git a/doc/src/Statistics/datafiles/automersenne.dat b/doc/src/Statistics/datafiles/automersenne.dat new file mode 100644 index 000000000..e60f67f00 --- /dev/null +++ b/doc/src/Statistics/datafiles/automersenne.dat @@ -0,0 +1,1000 @@ + 0 1.0000000 + 1 -0.00039850883 + 2 0.0075297237 + 3 -0.076614584 + 4 -0.011132200 + 5 -0.030989810 + 6 0.041812361 + 7 0.022850567 + 8 0.020375660 + 9 0.021512042 + 10 -0.038233298 + 11 -0.038677806 + 12 0.039211752 + 13 -0.011095491 + 14 -0.0044564329 + 15 -0.035488512 + 16 -0.036152110 + 17 -0.013012335 + 18 0.0020817233 + 19 -0.018777970 + 20 0.0015342904 + 21 -0.034707469 + 22 -0.0045328442 + 23 -0.022961177 + 24 -0.017504691 + 25 0.025258106 + 26 0.038105998 + 27 0.027601135 + 28 -0.0042234787 + 29 -0.048024298 + 30 0.046726347 + 31 0.033528879 + 32 0.010739866 + 33 0.048112058 + 34 0.074342844 + 35 -0.0035160703 + 36 -0.097559373 + 37 -0.0036536069 + 38 -0.021169986 + 39 -0.026994281 + 40 0.0011361472 + 41 0.00042746464 + 42 0.0094556642 + 43 0.035186489 + 44 0.012008881 + 45 -0.066250195 + 46 -0.0058885570 + 47 -0.0013470536 + 48 -0.020357559 + 49 0.013659362 + 50 -0.014194533 + 51 0.064474272 + 52 -0.049315732 + 53 0.0068481500 + 54 -0.0015570919 + 55 -0.041760952 + 56 0.047017059 + 57 0.0083207514 + 58 -0.018694355 + 59 -0.011955604 + 60 -0.018575782 + 61 -0.022540700 + 62 -0.0038630138 + 63 0.020907804 + 64 -0.024616327 + 65 0.11530859 + 66 -0.067525916 + 67 -0.0067679184 + 68 -0.032918810 + 69 -0.018539289 + 70 -0.084216861 + 71 0.019426308 + 72 -0.00074715715 + 73 0.037295939 + 74 0.029917292 + 75 -0.027646591 + 76 0.0044122468 + 77 0.060618328 + 78 0.051236343 + 79 -0.016390863 + 80 -0.0099063162 + 81 -0.00082734335 + 82 -0.014929417 + 83 -0.0091957452 + 84 0.026182610 + 85 0.0018857001 + 86 0.031396746 + 87 -0.045390441 + 88 -0.055549108 + 89 -0.052299358 + 90 -0.015109536 + 91 0.032190708 + 92 -0.026868033 + 93 0.021618018 + 94 0.0050115645 + 95 0.011128105 + 96 -0.055351651 + 97 0.023984066 + 98 0.025469199 + 99 0.024346935 + 100 -0.0096941685 + 101 -0.053075665 + 102 0.013409246 + 103 -0.0071408978 + 104 -0.012223020 + 105 -0.0016478848 + 106 0.021045263 + 107 0.046403350 + 108 0.0069554279 + 109 0.014708717 + 110 -0.00026539003 + 111 0.024845445 + 112 -0.0049231997 + 113 -0.0017220862 + 114 0.0045621071 + 115 -0.025192401 + 116 0.0028209063 + 117 -0.092084217 + 118 0.016684106 + 119 0.021835229 + 120 0.018734216 + 121 -0.036952674 + 122 -0.040815536 + 123 0.0015671699 + 124 0.010132854 + 125 0.046715210 + 126 -0.019535585 + 127 0.0044946246 + 128 0.026224408 + 129 -0.031652034 + 130 0.0094934804 + 131 -0.050688941 + 132 0.032252957 + 133 -0.0078409521 + 134 -0.011630815 + 135 -0.0087330534 + 136 0.039295136 + 137 0.045309325 + 138 -0.022338545 + 139 -0.021070208 + 140 0.039383224 + 141 -0.021284479 + 142 0.017074248 + 143 -5.1343887E-05 + 144 -0.045458822 + 145 -0.046122944 + 146 -0.035461130 + 147 -0.021600258 + 148 0.021437230 + 149 -0.0061934076 + 150 0.051751485 + 151 -0.025273260 + 152 -0.017574909 + 153 0.0067792040 + 154 0.0018960306 + 155 0.016680514 + 156 0.017963759 + 157 0.043943111 + 158 -0.0061690099 + 159 0.026941385 + 160 0.012216633 + 161 -0.019785526 + 162 0.033424289 + 163 0.0088921176 + 164 -0.043272488 + 165 -0.0073920701 + 166 -0.022735988 + 167 0.019831288 + 168 -0.022406716 + 169 -0.0034270288 + 170 0.038661913 + 171 0.0015083024 + 172 -0.022180629 + 173 -0.013418595 + 174 0.017036953 + 175 -0.012137219 + 176 0.0024605105 + 177 -0.033808969 + 178 -0.018628761 + 179 -0.027346843 + 180 -0.012037812 + 181 0.019175778 + 182 -0.066844592 + 183 0.040087724 + 184 0.0045306968 + 185 -0.0016694184 + 186 0.027033664 + 187 0.038560397 + 188 -0.0030428931 + 189 -0.014561361 + 190 0.058253381 + 191 -0.030927654 + 192 0.025350468 + 193 0.0081393072 + 194 -0.063504391 + 195 -0.048869374 + 196 0.061287388 + 197 0.0067389738 + 198 -0.014016816 + 199 0.017009642 + 200 0.013498226 + 201 0.066760012 + 202 -0.0078831629 + 203 -0.014284510 + 204 0.015891370 + 205 -0.0017101068 + 206 0.017561522 + 207 -0.042216875 + 208 0.0010105827 + 209 -0.039588987 + 210 -0.010623434 + 211 0.012428578 + 212 -0.044573760 + 213 0.080482020 + 214 0.036705368 + 215 -0.025105371 + 216 0.014592058 + 217 -0.017789069 + 218 0.038765246 + 219 -0.025885902 + 220 0.017043839 + 221 0.025500897 + 222 0.013596501 + 223 -0.019972053 + 224 -0.016509130 + 225 0.0027069680 + 226 0.014210155 + 227 0.017071603 + 228 0.0034474582 + 229 -0.019289487 + 230 0.032128214 + 231 -0.023774060 + 232 -0.0052932205 + 233 0.0043239495 + 234 0.049467003 + 235 -0.0019682175 + 236 -0.0057175850 + 237 -0.086458240 + 238 -0.010633932 + 239 0.049711797 + 240 -0.0051473949 + 241 -0.010502996 + 242 -0.059682111 + 243 -0.013875569 + 244 -0.050671337 + 245 -0.042970472 + 246 0.0032983110 + 247 0.022735402 + 248 0.051786772 + 249 -0.065820501 + 250 -0.014616146 + 251 -0.00053863927 + 252 0.036148981 + 253 -0.027105785 + 254 -0.0033587844 + 255 0.012587625 + 256 -0.032473465 + 257 0.0045156150 + 258 0.0025041393 + 259 0.022749381 + 260 -0.011942502 + 261 0.040564347 + 262 0.019153347 + 263 0.042605069 + 264 0.026921673 + 265 -0.017759274 + 266 0.014206232 + 267 -0.0069809192 + 268 -0.0011496859 + 269 0.0063889545 + 270 -0.0075013878 + 271 -0.046046299 + 272 -0.013068379 + 273 -0.0062615689 + 274 0.0056876531 + 275 -0.025717888 + 276 -0.0072959257 + 277 -0.025127295 + 278 0.044797986 + 279 -0.011806658 + 280 -0.011621526 + 281 0.016598971 + 282 -0.023535958 + 283 -0.0033071575 + 284 0.0039764024 + 285 0.032759776 + 286 -0.0072796509 + 287 -0.022481839 + 288 -0.038837342 + 289 -0.035360598 + 290 0.047625475 + 291 0.0029441419 + 292 0.015363520 + 293 -0.032890305 + 294 0.00014385169 + 295 0.0093506211 + 296 0.042240558 + 297 0.016496134 + 298 -0.0046955917 + 299 0.026925193 + 300 -0.047930589 + 301 0.022580448 + 302 -0.041639190 + 303 0.031523493 + 304 -0.029478734 + 305 0.021526621 + 306 -0.022861985 + 307 0.0016559278 + 308 0.0028637402 + 309 0.026062242 + 310 0.0015438354 + 311 0.035893343 + 312 0.025777036 + 313 0.020541070 + 314 -0.056851732 + 315 0.034370506 + 316 -0.034130810 + 317 0.0032066005 + 318 -0.028514122 + 319 -0.016266227 + 320 0.014115845 + 321 0.0098882228 + 322 -0.046938829 + 323 0.030429223 + 324 0.049126829 + 325 -0.020672407 + 326 -0.0066021677 + 327 0.0017879728 + 328 -0.026402363 + 329 -0.0040738627 + 330 -0.025430820 + 331 -0.027299590 + 332 -0.020525149 + 333 0.040032531 + 334 0.0041808025 + 335 -0.0036519573 + 336 -0.0027050634 + 337 -0.024097320 + 338 0.019582987 + 339 0.045558977 + 340 0.0067800239 + 341 -0.0042827926 + 342 -0.011481723 + 343 0.041870674 + 344 -0.0078637289 + 345 0.0047561047 + 346 0.042687334 + 347 -0.0053904870 + 348 0.0066880645 + 349 0.00076997894 + 350 0.015016180 + 351 -0.046314612 + 352 -0.0035992390 + 353 0.018856659 + 354 0.029241951 + 355 -0.023837018 + 356 -0.0092689546 + 357 -0.027062945 + 358 0.017439465 + 359 -0.011280277 + 360 -0.020408789 + 361 0.029352588 + 362 -0.0065214405 + 363 0.0019510915 + 364 0.030417934 + 365 -0.021614425 + 366 0.040519121 + 367 -0.014387146 + 368 0.022631461 + 369 -0.037162875 + 370 0.057832805 + 371 -0.0012082448 + 372 -0.025176406 + 373 0.021927684 + 374 0.036809227 + 375 0.011755410 + 376 0.039719074 + 377 -0.0057936908 + 378 -0.014894587 + 379 -0.030001782 + 380 -0.0068359635 + 381 -0.027206809 + 382 0.0033147001 + 383 0.019952090 + 384 -0.016189684 + 385 -0.014190482 + 386 0.013660921 + 387 0.020810808 + 388 0.036953804 + 389 0.0068815608 + 390 -0.0051748463 + 391 -0.035329903 + 392 0.017989289 + 393 0.0017726599 + 394 -0.028362238 + 395 0.0076066270 + 396 0.012315135 + 397 0.019602391 + 398 0.0054637460 + 399 0.010180415 + 400 0.010715276 + 401 0.032347243 + 402 -0.011310452 + 403 -0.019241395 + 404 0.046814252 + 405 0.0033537765 + 406 -0.0080705101 + 407 -0.034532463 + 408 -0.013596282 + 409 -0.030167978 + 410 0.021618277 + 411 0.023256360 + 412 -0.024011435 + 413 0.0099447977 + 414 0.0051984929 + 415 0.0032058618 + 416 -0.018386174 + 417 0.091155277 + 418 0.0038552629 + 419 -0.011819495 + 420 -0.019909822 + 421 0.0049708514 + 422 0.030180281 + 423 0.019513707 + 424 0.010365884 + 425 -0.021038696 + 426 0.012563095 + 427 -0.0078607670 + 428 -0.034763368 + 429 -0.013946109 + 430 0.0018269461 + 431 0.034874965 + 432 -0.027224239 + 433 0.015590005 + 434 -0.019383632 + 435 0.0038593809 + 436 -0.0056616145 + 437 -0.020352622 + 438 0.013123802 + 439 -0.025004282 + 440 0.0018341696 + 441 -0.039770439 + 442 0.00087601823 + 443 -0.0066024767 + 444 -0.0097114643 + 445 -0.011356804 + 446 -0.0071475429 + 447 0.035558002 + 448 0.021950580 + 449 0.0071515180 + 450 -0.019641713 + 451 -0.012505286 + 452 0.00058888713 + 453 -0.0032990818 + 454 0.0055603481 + 455 -0.012210133 + 456 -0.036752880 + 457 -0.0010291934 + 458 0.013480362 + 459 0.0011063715 + 460 0.015099980 + 461 0.010381880 + 462 -0.012656722 + 463 -0.034560636 + 464 0.046300125 + 465 -0.012855784 + 466 -0.0069566210 + 467 -0.014963271 + 468 0.0071160662 + 469 -0.023183606 + 470 -0.0054130250 + 471 -0.029441744 + 472 -0.013208213 + 473 0.0058226359 + 474 0.0073592316 + 475 -0.026501589 + 476 0.0086787125 + 477 0.011144501 + 478 -0.017992368 + 479 -0.020137351 + 480 -0.00050191872 + 481 0.032606384 + 482 0.028294831 + 483 0.0069425809 + 484 -0.0025918567 + 485 -0.030876138 + 486 0.028514704 + 487 -0.030215330 + 488 -0.022505299 + 489 -0.0045540661 + 490 -0.021775574 + 491 0.00038863032 + 492 -0.019469717 + 493 -0.026643162 + 494 0.015456519 + 495 0.0060890381 + 496 0.0072811940 + 497 -0.020604861 + 498 0.0058975829 + 499 0.0077669157 + 500 -0.014203543 + 501 0.027618735 + 502 0.010987245 + 503 0.016146396 + 504 -0.046623086 + 505 -0.021404063 + 506 -0.0071026109 + 507 0.014316482 + 508 0.050116141 + 509 -0.016053678 + 510 0.00037589453 + 511 0.024107855 + 512 -0.020410421 + 513 -0.0039090131 + 514 0.018030706 + 515 -0.019683084 + 516 0.020915648 + 517 -0.012283794 + 518 -0.072770967 + 519 -0.0020466698 + 520 0.033620334 + 521 0.0026312146 + 522 -0.019340110 + 523 -0.010093083 + 524 0.032372824 + 525 0.023358211 + 526 0.019563085 + 527 -0.023193643 + 528 -0.023912984 + 529 0.028365424 + 530 -0.017833524 + 531 0.049366485 + 532 0.0033034666 + 533 0.036273527 + 534 -0.034107635 + 535 0.027272035 + 536 -0.012223060 + 537 0.047632726 + 538 -0.010089959 + 539 -0.0028197945 + 540 0.00065234228 + 541 0.0032028460 + 542 -0.014032424 + 543 -0.0038940314 + 544 -0.014471639 + 545 0.010245881 + 546 -0.031846228 + 547 -0.018379583 + 548 -0.012049365 + 549 -0.010966222 + 550 -0.013878160 + 551 0.016589271 + 552 -0.021916867 + 553 0.0087858521 + 554 0.027621607 + 555 -0.018883327 + 556 0.0055198510 + 557 -0.0060233484 + 558 0.0065385887 + 559 0.049948053 + 560 0.022709131 + 561 0.0083870584 + 562 0.029122276 + 563 0.030653719 + 564 -0.0038621603 + 565 -0.0077334777 + 566 0.0013917389 + 567 0.025550810 + 568 -0.011855303 + 569 -0.030284018 + 570 -0.0073176655 + 571 0.017890083 + 572 0.0023946457 + 573 0.0082304715 + 574 -0.014303315 + 575 0.020850567 + 576 -0.010499909 + 577 -0.028686313 + 578 0.00029063430 + 579 0.0058666114 + 580 -0.021785486 + 581 -0.0096116460 + 582 -0.010938017 + 583 -0.028888331 + 584 -0.00055043544 + 585 0.020775413 + 586 0.019325640 + 587 0.021697931 + 588 -0.031392595 + 589 -0.033268006 + 590 -0.039920469 + 591 -0.015858921 + 592 -0.0010340578 + 593 0.028674925 + 594 0.027232115 + 595 -0.023133588 + 596 -0.031392641 + 597 -0.0033751189 + 598 -0.043930110 + 599 -0.022065423 + 600 0.018575533 + 601 0.034786503 + 602 0.022188960 + 603 -0.018722519 + 604 -0.0033538481 + 605 -0.0019135364 + 606 0.034412541 + 607 0.0063933819 + 608 0.0043941754 + 609 0.0031151635 + 610 0.0059123942 + 611 -0.022565534 + 612 -0.00012048878 + 613 0.018358466 + 614 0.0058849744 + 615 -0.010446412 + 616 0.027175242 + 617 0.015094912 + 618 0.045011644 + 619 0.0051040681 + 620 -0.015088009 + 621 0.00060846204 + 622 -0.023562337 + 623 -0.014687213 + 624 -0.0019508561 + 625 -0.027180327 + 626 0.0053059608 + 627 0.026911053 + 628 0.019832294 + 629 -0.014316090 + 630 -0.013579816 + 631 -0.0074524339 + 632 0.0023041549 + 633 -0.0048902581 + 634 0.014748700 + 635 0.015828956 + 636 0.011651372 + 637 -0.028449219 + 638 0.018893515 + 639 0.022430458 + 640 0.0020572850 + 641 -0.047852942 + 642 -0.030971596 + 643 -0.011280746 + 644 0.035021969 + 645 -0.011408873 + 646 -0.024270635 + 647 -0.024445963 + 648 -0.0080845450 + 649 0.020187707 + 650 0.010225120 + 651 0.026541612 + 652 -0.0094225257 + 653 -0.028655708 + 654 -0.028771259 + 655 -0.054234804 + 656 0.027811754 + 657 0.012886491 + 658 0.027248276 + 659 0.0014272864 + 660 -0.0038342847 + 661 0.012264597 + 662 0.0078784196 + 663 -0.014772008 + 664 -0.0021303419 + 665 -0.0023548533 + 666 -0.010001963 + 667 0.010265614 + 668 -0.016169868 + 669 0.00068463391 + 670 0.0019238310 + 671 0.00054326870 + 672 -0.0016790464 + 673 0.020675073 + 674 -0.030547720 + 675 -0.016952851 + 676 -0.0043647773 + 677 0.0089803777 + 678 -0.0011795400 + 679 -0.0026837760 + 680 -0.0066145568 + 681 0.016555175 + 682 0.0045479616 + 683 0.021909635 + 684 -0.0070725479 + 685 0.011858958 + 686 0.0063636282 + 687 0.013213975 + 688 0.0086205024 + 689 0.0034084325 + 690 -0.012270779 + 691 0.0018423414 + 692 0.013453993 + 693 0.0050505880 + 694 0.022730152 + 695 0.010188657 + 696 -0.0051335682 + 697 -0.0020429929 + 698 -0.0046473046 + 699 -0.0064359202 + 700 0.027437169 + 701 0.0021062070 + 702 -0.0092277785 + 703 -0.021793309 + 704 -0.0093357745 + 705 0.0096844781 + 706 -0.017812423 + 707 0.019216536 + 708 -0.011751425 + 709 0.031306164 + 710 -0.031003733 + 711 0.0045908435 + 712 0.022197245 + 713 0.044356823 + 714 -0.011279340 + 715 -0.00063210160 + 716 0.022619262 + 717 -0.0048696777 + 718 2.7710142E-05 + 719 0.0093208258 + 720 0.011151741 + 721 0.042734701 + 722 0.00082345063 + 723 -0.0064815085 + 724 -0.011162417 + 725 0.019371743 + 726 -0.0072227478 + 727 0.036856608 + 728 -0.0046597009 + 729 0.0093292407 + 730 -0.029881720 + 731 -0.022544055 + 732 0.0067547008 + 733 -0.012852388 + 734 0.017629149 + 735 0.022035951 + 736 0.0039029374 + 737 0.0038631069 + 738 -0.0064078979 + 739 -0.025045521 + 740 -0.030156360 + 741 0.0022182783 + 742 0.0056079074 + 743 -0.0057498694 + 744 -0.0038137357 + 745 -0.0019111597 + 746 0.033244546 + 747 0.0039995162 + 748 -0.028431035 + 749 -0.0084654091 + 750 -0.017920230 + 751 0.0091895467 + 752 -0.047359823 + 753 0.0053614276 + 754 0.0045515416 + 755 0.0063096922 + 756 -0.0051722096 + 757 -0.0090646204 + 758 -0.015375531 + 759 0.010469262 + 760 0.020619059 + 761 0.019446803 + 762 0.015360235 + 763 0.026556079 + 764 -0.014850757 + 765 0.00043243229 + 766 -0.019207935 + 767 0.0091227560 + 768 -0.0039674501 + 769 0.012677223 + 770 0.0059764792 + 771 -0.0019081743 + 772 0.00052190352 + 773 -0.023890356 + 774 -0.0093265644 + 775 -0.0027615351 + 776 -0.010170966 + 777 -0.026209697 + 778 -0.0048646684 + 779 0.0030630937 + 780 -0.0039676134 + 781 -0.0088915791 + 782 -0.014279751 + 783 -0.019704165 + 784 -0.019209887 + 785 0.0061373699 + 786 0.020227079 + 787 0.028729569 + 788 0.012112855 + 789 0.017006455 + 790 -0.0033648302 + 791 0.00030061169 + 792 0.0012253500 + 793 0.023771468 + 794 0.0046375357 + 795 -0.011610943 + 796 -0.018452313 + 797 0.0083301500 + 798 0.0037865728 + 799 0.021777084 + 800 0.015769708 + 801 -0.014965320 + 802 -0.0017419876 + 803 -0.033476334 + 804 -0.030830255 + 805 -0.021870028 + 806 0.0074449924 + 807 0.0013585510 + 808 -0.026724716 + 809 -0.0089502027 + 810 0.0047979763 + 811 -0.013332303 + 812 -0.0064698535 + 813 -0.012448045 + 814 -0.0047413069 + 815 -0.0038978876 + 816 0.016833818 + 817 0.014603311 + 818 0.020443440 + 819 0.028826407 + 820 0.013123271 + 821 0.0091723388 + 822 -0.011255640 + 823 -0.0012867583 + 824 0.0074609431 + 825 0.0079598103 + 826 0.00055143813 + 827 0.0073569375 + 828 0.027685980 + 829 -0.0079481602 + 830 -0.015245658 + 831 -0.016568012 + 832 0.014166210 + 833 -0.017490318 + 834 0.0068715750 + 835 -0.010515617 + 836 -0.023031012 + 837 -0.0039414950 + 838 -0.021860724 + 839 0.00029086251 + 840 0.025457721 + 841 0.018815201 + 842 -0.014766348 + 843 -0.020393215 + 844 0.025794784 + 845 -0.0082851118 + 846 0.00013439103 + 847 0.0045753418 + 848 0.0080142804 + 849 -0.0056724468 + 850 -0.0053680915 + 851 -5.0186689E-05 + 852 -0.0032442615 + 853 0.0060424922 + 854 -0.017920096 + 855 -0.018203065 + 856 -0.0012571564 + 857 0.0078961550 + 858 -0.014874879 + 859 -0.0071654637 + 860 -0.012561047 + 861 -0.022924028 + 862 -0.0010177375 + 863 -0.016723859 + 864 -0.0041250596 + 865 -0.0099232278 + 866 -0.012007622 + 867 -0.010921927 + 868 -0.010611169 + 869 -0.0064454145 + 870 -0.012373089 + 871 0.016952049 + 872 0.0052607349 + 873 -0.012175133 + 874 0.0054016304 + 875 0.0092675688 + 876 -0.020126948 + 877 -0.015452958 + 878 0.0042092441 + 879 0.011057301 + 880 -0.014070768 + 881 0.0039739742 + 882 -0.010045862 + 883 -0.0024596864 + 884 -0.0025652609 + 885 -0.0046209025 + 886 0.011785118 + 887 -0.0069725970 + 888 -0.017379060 + 889 -0.024952371 + 890 -0.0037120511 + 891 0.0062956522 + 892 -0.015428430 + 893 0.014351959 + 894 -0.0074217601 + 895 -0.00016310749 + 896 0.00035103207 + 897 0.029621595 + 898 0.0095742916 + 899 0.012108004 + 900 -0.0013415321 + 901 0.0040298110 + 902 0.0081285013 + 903 0.0010463687 + 904 -0.013433386 + 905 0.010687478 + 906 -0.0016356944 + 907 -0.012058460 + 908 -0.0036622641 + 909 0.0021945389 + 910 -0.011631748 + 911 -0.012051364 + 912 0.013892661 + 913 -4.9505611E-05 + 914 -0.0016314468 + 915 7.4078423E-05 + 916 -0.011434937 + 917 -0.0069514616 + 918 0.016964461 + 919 0.012711562 + 920 -0.011803641 + 921 0.0059246183 + 922 0.0084047975 + 923 0.00085369214 + 924 0.013476117 + 925 0.017477181 + 926 0.013697700 + 927 0.0078041173 + 928 0.0081624940 + 929 0.0079026863 + 930 0.018032525 + 931 -0.0021333952 + 932 -0.0089801782 + 933 -0.0043425879 + 934 -0.0012769913 + 935 0.0031995993 + 936 -0.0085301623 + 937 0.0016565121 + 938 -0.016048639 + 939 0.0041939700 + 940 -0.0042789437 + 941 0.0040290875 + 942 -0.0078984469 + 943 -0.0031005041 + 944 -0.010489387 + 945 -0.014048278 + 946 0.0080731664 + 947 0.0072791508 + 948 0.017732252 + 949 0.0070456969 + 950 0.0048539841 + 951 -0.0063464819 + 952 0.0071559674 + 953 -0.00076997791 + 954 -0.0060856751 + 955 0.0052101682 + 956 0.0076901419 + 957 -0.0027116074 + 958 -0.0015772656 + 959 0.011581456 + 960 0.0052794831 + 961 -0.0076978630 + 962 -0.0081319124 + 963 -0.0092991557 + 964 -0.012073101 + 965 -0.013250383 + 966 -0.0036785721 + 967 0.0050898879 + 968 0.0028822611 + 969 -0.010407446 + 970 0.0011417214 + 971 -0.0074299748 + 972 -0.0015352307 + 973 -0.0057842566 + 974 0.011529249 + 975 0.0087991237 + 976 0.00030885135 + 977 0.0036026354 + 978 -0.0014831089 + 979 0.00037216609 + 980 -0.0014774795 + 981 0.0074766204 + 982 0.0072999330 + 983 0.0055104880 + 984 0.0020526707 + 985 -0.0042872281 + 986 -0.0023896493 + 987 -0.0060384437 + 988 -0.0070743121 + 989 -0.0062837043 + 990 -0.0039570614 + 991 -0.0034175474 + 992 0.0015737944 + 993 -6.8705654E-05 + 994 0.0059629554 + 995 0.0020717095 + 996 0.0023956585 + 997 -0.00034251251 + 998 0.0015951958 + 999 0.00031542955 diff --git a/doc/src/Statistics/datafiles/plot.py b/doc/src/Statistics/datafiles/plot.py new file mode 100644 index 000000000..ecea29a77 --- /dev/null +++ b/doc/src/Statistics/datafiles/plot.py @@ -0,0 +1,16 @@ +import numpy as np +from matplotlib import pyplot as plt +# Load in data file +data = np.loadtxt("autocor.dat") +data1 = np.loadtxt("automersenne.dat") +# Make arrays containing x-axis and binding energies as function of A +x = data[:,0] +corr = data[:,1] +corr2 = data1[:,1] +plt.plot(x, corr ,'ro', x, corr2, 'b') +plt.axis([0,1000,-0.2, 1.1]) +plt.xlabel(r'$d$') +plt.ylabel(r'$C_d$') +plt.title(r'autocorrelation function for RNG') +plt.savefig('autocorr.pdf') +plt.show() diff --git a/doc/src/Statistics/make.sh b/doc/src/Statistics/make.sh new file mode 100755 index 000000000..edfed7c87 --- /dev/null +++ b/doc/src/Statistics/make.sh @@ -0,0 +1,118 @@ +#!/bin/sh +set -x + +function system { + "$@" + if [ $? -ne 0 ]; then + echo "make.sh: unsuccessful command $@" + echo "abort!" + exit 1 + fi +} + +if [ $# -eq 0 ]; then +echo 'bash make.sh slides1|slides2' +exit 1 +fi + +name=$1 +rm -f *.tar.gz + +opt="--encoding=utf-8" +# Note: Makefile examples contain constructions like ${PROG} which +# looks like Mako constructions, but they are not. Use --no_mako +# to turn off Mako processing. +opt="--no_mako" + +rm -f *.aux + + +html=${name}-reveal +system doconce format html $name --pygments_html_style=perldoc --keep_pygments_html_bg --html_links_in_new_window --html_output=$html $opt +system doconce slides_html $html reveal --html_slide_theme=beige + +# Plain HTML documents + +html=${name}-solarized +system doconce format html $name --pygments_html_style=perldoc --html_style=solarized3 --html_links_in_new_window --html_output=$html $opt +system doconce split_html $html.html --method=space10 + +html=${name} +system doconce format html $name --pygments_html_style=default --html_style=bloodish --html_links_in_new_window --html_output=$html $opt +system doconce split_html $html.html --method=space10 + +# Bootstrap style +html=${name}-bs +system doconce format html $name --html_style=bootstrap --pygments_html_style=default --html_admon=bootstrap_panel --html_output=$html $opt +system doconce split_html $html.html --method=split --pagination --nav_button=bottom + +# IPython notebook +system doconce format ipynb $name $opt + +# LaTeX Beamer slides +beamertheme=red_plain +system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt +system doconce ptex2tex $name envir=minted +# Add special packages +doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex +system doconce slides_beamer $name --beamer_slide_theme=$beamertheme +system pdflatex -shell-escape ${name} +system pdflatex -shell-escape ${name} +cp $name.pdf ${name}-beamer.pdf +cp $name.tex ${name}-beamer.tex + +# Handouts +system doconce format pdflatex $name --latex_title_layout=beamer --latex_table_format=footnotesize $opt +system doconce ptex2tex $name envir=minted +# Add special packages +doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex +system doconce slides_beamer $name --beamer_slide_theme=red_shadow --handout +system pdflatex -shell-escape $name +pdflatex -shell-escape $name +pdflatex -shell-escape $name +pdfnup --nup 2x3 --frame true --delta "1cm 1cm" --scale 0.9 --outfile ${name}-beamer-handouts2x3.pdf ${name}.pdf +rm -f ${name}.pdf + +# Ordinary plain LaTeX document +rm -f *.aux # important after beamer +system doconce format pdflatex $name --minted_latex_style=trac --latex_admon=paragraph $opt +system doconce ptex2tex $name envir=minted +# Add special packages +doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex +doconce replace 'section{' 'section*{' $name.tex +pdflatex -shell-escape $name +pdflatex -shell-escape $name +mv -f $name.pdf ${name}-minted.pdf +cp $name.tex ${name}-plain-minted.tex + + + +# Publish +dest=../../pub +if [ ! -d $dest/$name ]; then +mkdir $dest/$name +mkdir $dest/$name/pdf +mkdir $dest/$name/html +mkdir $dest/$name/ipynb +fi +cp ${name}*.pdf $dest/$name/pdf +cp -r ${name}*.html ._${name}*.html reveal.js $dest/$name/html + +# Figures: cannot just copy link, need to physically copy the files +if [ -d fig-${name} ]; then +if [ ! -d $dest/$name/html/fig-$name ]; then +mkdir $dest/$name/html/fig-$name +fi +cp -r fig-${name}/* $dest/$name/html/fig-$name +fi + +cp ${name}.ipynb $dest/$name/ipynb +ipynb_tarfile=ipynb-${name}-src.tar.gz +if [ ! -f ${ipynb_tarfile} ]; then +cat > README.txt <