From bc1dc849df925284a38d23e33433e08fbdbf8b49 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Tue, 25 Aug 2020 10:53:39 +0200 Subject: [PATCH] test --- doc/src/Statistics/Statistics.do.txt~ | 1867 ------------------------- 1 file changed, 1867 deletions(-) delete mode 100644 doc/src/Statistics/Statistics.do.txt~ diff --git a/doc/src/Statistics/Statistics.do.txt~ b/doc/src/Statistics/Statistics.do.txt~ deleted file mode 100644 index f282d896e..000000000 --- a/doc/src/Statistics/Statistics.do.txt~ +++ /dev/null @@ -1,1867 +0,0 @@ -TITLE: Data Analysis and Machine Learning: Elements of Probability Theory and Statistical Data Analysis -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 -===== Domains and probabilities ===== -!bblock -Consider the following simple example, namely the tossing of two 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 -===== Tossing the dice ===== -!bblock -The numbers in the domain are the outcomes of the physical process of tossing say two 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 probability distribution function(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 -For $a=0$ and $b=1$ we have -!bt -\[ -\begin{array}{ll} -p(x)dx = dx & \in [0,1]. -\end{array} -\] -!et -The latter distribution is used to generate random numbers. For other PDFs, one needs normally a mapping from this distribution to say for example the exponential distribution. -!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 \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 - 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 -===== 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=\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 -===== Covariance ===== -!bc pycod -# Importing various packages -from math import exp, sqrt -from random import random, seed -import numpy as np -import matplotlib.pyplot as plt - -def covariance(x, y, n): - sum = 0.0 - mean_x = np.mean(x) - mean_y = np.mean(y) - for i in range(0, n): - sum += (x[(i)]-mean_x)*(y[i]-mean_y) - return sum/n - -n = 10 - -x=np.random.normal(size=n) -y = 4+3*x+np.random.normal(size=n) -covxy = covariance(x,y,n) -print(covxy) -z = np.vstack((x, y)) -c = np.cov(z.T) - -print(c) - -!ec - - - -!split -===== Meet the covariance, uncorrelated events ===== -!bblock - -Consider the stochastic variables $X_i$ and $X_j$, ($i\neq j$). We have -!bt -\begin{align*} -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 (assuming $i \neq j$), we have that -!bt -\[ -\langle x_i x_j\rangle = \langle x_i\rangle\langle x_j\rangle, -\] -!et -leading to -!bt -\[ -Cov(X_i, X_j) = 0 \hspace{0.1cm} (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 -======= Which RNG should I use? ======= -!bblock -* 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}$. -* Add RNGs in Python - -!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 -===== Why blocking? ===== -!bblock Statistical analysis - * Monte Carlo simulations can be treated as *computer experiments* - * The results can be analysed with the same statistical tools as we would use analysing experimental data. - * As in all experiments, we are looking for expectation values and an estimate of how accurate they are, i.e., possible sources for errors. - -A very good article which explains blocking is H. Flyvbjerg and H. G. Petersen, *Error estimates on averages of correlated data*, "Journal of Chemical Physics 91, 461-466 (1989)":"http://scitation.aip.org/content/aip/journal/jcp/91/1/10.1063/1.457480". - -!eblock - - - -!split -===== Why blocking? ===== -!bblock Statistical analysis - * As in other experiments, Monte Carlo experiments have two classes of errors: - * Statistical errors - * Systematical errors - * Statistical errors can be estimated using standard tools from statistics - * Systematical errors are method specific and must be treated differently from case to case. (In VMC a common source is the step length or time step in importance sampling) -!eblock - - -!split -===== Code to demonstrate the calculation of the autocorrelation function ===== -The following code computes the autocorrelation function, the covariance and the standard deviation -for standard RNG. -The "following file":"https://github.com/CompPhysics/ComputationalPhysics2/tree/gh-pages/doc/Programs/LecturePrograms/programs/Blocking/autocorrelation.cpp" gives the code. -!bc cppcod -// This function computes the autocorrelation function for -// the Mersenne random number generator with a uniform distribution -#include -#include -#include -#include -#include -#include -#include -#include -using namespace std; -using namespace arma; -// output file -ofstream ofile; - -// Main function begins here -int main(int argc, char* argv[]) -{ - int MonteCarloCycles; - string filename; - if (argc > 1) { - filename=argv[1]; - MonteCarloCycles = atoi(argv[2]); - string fileout = filename; - string argument = to_string(MonteCarloCycles); - fileout.append(argument); - ofile.open(fileout); - } - - // 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 - // compute the covariance and the correlation function - - vec X = zeros(MonteCarloCycles); - double MCint = 0.; double MCintsqr2=0.; - 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); - for (int i = 0; i < MonteCarloCycles; i++){ - double x = RandomNumberGenerator(gen); - X(i) = x; - MCint += x; - MCintsqr2 += x*x; - } - double Mean = MCint/((double) MonteCarloCycles ); - MCintsqr2 = MCintsqr2/((double) MonteCarloCycles ); - double STDev = sqrt(MCintsqr2-Mean*Mean); - double Variance = MCintsqr2-Mean*Mean; - // Write mean value and variance - cout << " Sample variance= " << Variance << " Mean value = " << Mean << endl; - // Now we compute the autocorrelation function - vec autocorrelation = zeros(MonteCarloCycles); - for (int j = 0; j < MonteCarloCycles; j++){ - double sum = 0.0; - for (int k = 0; k < (MonteCarloCycles-j); k++){ - sum += (X(k)-Mean)*(X(k+j)-Mean); - } - autocorrelation(j) = sum/Variance/((double) MonteCarloCycles ); - ofile << setiosflags(ios::showpoint | ios::uppercase); - ofile << setw(15) << setprecision(8) << j; - ofile << setw(15) << setprecision(8) << autocorrelation(j) << endl; - } - // Now compute the exact covariance using the autocorrelation function - double Covariance = 0.0; - for (int j = 0; j < MonteCarloCycles; j++){ - Covariance += autocorrelation(j); - } - Covariance *= 2.0/((double) MonteCarloCycles); - // Compute now the total variance, including the covariance, and obtain the standard deviation - double TotalVariance = (Variance/((double) MonteCarloCycles ))+Covariance; - cout << "Covariance =" << Covariance << "Totalvariance= " << TotalVariance << "Sample Variance/n= " << (Variance/((double) MonteCarloCycles )) << endl; - cout << " STD from sample variance= " << sqrt(Variance/((double) MonteCarloCycles )) << " STD with covariance = " << sqrt(TotalVariance) << endl; - - ofile.close(); // close output file - return 0; -} // end of main program - - -!ec - - -!split -===== What is blocking? ===== -!bblock Blocking - * Say that we have a set of samples from a Monte Carlo experiment - * Assuming (wrongly) that our samples are uncorrelated our best estimate of the standard deviation of the mean $\langle \mathbf{M}\rangle$ is given by -!bt -\[ -\sigma=\sqrt{\frac{1}{n}\left(\langle \mathbf{M}^2\rangle-\langle \mathbf{M}\rangle^2\right)} -\] -!et - * If the samples are correlated we can rewrite our results to show that -!bt -\[ -\sigma=\sqrt{\frac{1+2\tau/\Delta t}{n}\left(\langle \mathbf{M}^2\rangle-\langle \mathbf{M}\rangle^2\right)} -\] -!et - where $\tau$ is the correlation time (the time between a sample and the next uncorrelated sample) and $\Delta t$ is time between each sample -!eblock - -!split -===== What is blocking? ===== -!bblock Blocking - * If $\Delta t\gg\tau$ our first estimate of $\sigma$ still holds - * Much more common that $\Delta t<\tau$ - * In the method of data blocking we divide the sequence of samples into blocks - * We then take the mean $\langle \mathbf{M}_i\rangle$ of block $i=1\ldots n_{blocks}$ to calculate the total mean and variance - * The size of each block must be so large that sample $j$ of block $i$ is not correlated with sample $j$ of block $i+1$ - * The correlation time $\tau$ would be a good choice -!eblock - -!split -===== What is blocking? ===== -!bblock Blocking - * Problem: We don't know $\tau$ or it is too expensive to compute - * Solution: Make a plot of std. dev. as a function of blocksize - * The estimate of std. dev. of correlated data is too low $\to$ the error will increase with increasing block size until the blocks are uncorrelated, where we reach a plateau - * When the std. dev. stops increasing the blocks are uncorrelated -!eblock - -!split -===== Implementation ===== -!bblock - * Do a Monte Carlo simulation, storing all samples to file - * Do the statistical analysis on this file, independently of your Monte Carlo program - * Read the file into an array - * Loop over various block sizes - * For each block size $n_b$, loop over the array in steps of $n_b$ taking the mean of elements $i n_b,\ldots,(i+1) n_b$ - * Take the mean and variance of the resulting array - * Write the results for each block size to file for later - analysis -!eblock - - - - - -!split -===== Actual implementation with code, main function ===== -When the file gets large, it can be useful to write your data in binary mode instead of ascii characters. -The "following python file":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py" reads data from file with the output from every Monte Carlo cycle. -!bc pycod -# Blocking - @timeFunction - def blocking(self, blockSizeMax = 500): - blockSizeMin = 1 - - self.blockSizes = [] - self.meanVec = [] - self.varVec = [] - - for i in range(blockSizeMin, blockSizeMax): - if(len(self.data) % i != 0): - pass#continue - blockSize = i - meanTempVec = [] - varTempVec = [] - startPoint = 0 - endPoint = blockSize - - while endPoint <= len(self.data): - meanTempVec.append(np.average(self.data[startPoint:endPoint])) - startPoint = endPoint - endPoint += blockSize - mean, var = np.average(meanTempVec), np.var(meanTempVec)/len(meanTempVec) - self.meanVec.append(mean) - self.varVec.append(var) - self.blockSizes.append(blockSize) - - self.blockingAvg = np.average(self.meanVec[-200:]) - self.blockingVar = (np.average(self.varVec[-200:])) - self.blockingStd = np.sqrt(self.blockingVar) - -!ec - - - - -!split -===== The Bootstrap method ===== - -The Bootstrap resampling method is also very popular. It is very simple: - -o Start with your sample of measurements and compute the sample variance and the mean values -o Then start again but pick in a random way the numbers in the sample and recalculate the mean and the sample variance. -o Repeat this $K$ times. - -It can be shown, see the article by "Efron":"https://projecteuclid.org/download/pdf_1/euclid.aos/1176344552" -that it produces the correct standard deviation. - -This method is very useful for small ensembles of data points. - -!split -===== Bootstrapping ===== -Given a set of $N$ data, assume that we are interested in some -observable $\theta$ which may be estimated from that set. This observable can also be for example the result of a fit based on all $N$ raw data. -Let us call the value of the observable obtained from the original -data set $\hat{\theta}$. One recreates from the sample repeatedly -other samples by choosing randomly $N$ data out of the original set. -This costs essentially nothing, since we just recycle the original data set for the building of new sets. - -!split -===== Bootstrapping, recipe ===== -Let us assume we have done this $K$ times and thus have $K$ sets of $N$ -data values each. -Of course some values will enter more than once in the new sets. For each of these sets one computes the observable $\theta$ resulting in values $\theta_k$ with $k = 1,...,K$. Then one determines -!bt -\[ -\tilde{\theta} = \frac{1}{K} \sum_{k=1}^K \theta_k, -\] -!et -and -!bt -\[ -sigma^2_{\tilde{\theta}} = \frac{1}{K} \sum_{k=1}^K \left(\theta_k-\tilde{\theta}\right)^2. -\] -!et - -These are estimators for $\angle\theta\rangle$ and its variance. They are not unbiased and therefore -$\tilde{\theta}\neq\hat{\theta}$ for finite K. - -The difference is called bias and gives an idea on how far away the result may be from -the true $\angle\theta\rangle$. As final result for the observable one quotes $\angle\theta\rangle = \tilde{\theta} \pm \sigma_{\tilde{\theta}}$ . - - -!split -===== Bootstrapping, "code":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py" ===== -!bc -# Bootstrap - @timeFunction - def bootstrap(self, nBoots = 1000): - bootVec = np.zeros(nBoots) - for k in range(0,nBoots): - bootVec[k] = np.average(np.random.choice(self.data, len(self.data))) - self.bootAvg = np.average(bootVec) - self.bootVar = np.var(bootVec) - self.bootStd = np.std(bootVec) -!ec - -!split -===== Jackknife, "code":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py" ===== -!bc -# Jackknife - @timeFunction - def jackknife(self): - jackknVec = np.zeros(len(self.data)) - for k in range(0,len(self.data)): - jackknVec[k] = np.average(np.delete(self.data, k)) - self.jackknAvg = self.avg - (len(self.data) - 1) * (np.average(jackknVec) - self.avg) - self.jackknVar = float(len(self.data) - 1) * np.var(jackknVec) - self.jackknStd = np.sqrt(self.jackknVar) -!ec - -