For a linear fit we don't need to invert a matrix!! +Defining $$ -\gamma = \sum_{i=0}^{n-1}\frac{n-1}{\sigma_i^2}, +\gamma = \sum_{i=0}^{n-1}\frac{1}{\sigma_i^2}, $$ @@ -231,7 +232,7 @@ $$ \gamma_{xy} = \sum_{i=0}^{n-1}\frac{y_ix_{i}}{\sigma_i^2}, $$ -and show that +we obtain $$ \beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, $$ @@ -241,7 +242,7 @@ $$ $$
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. +This approach (different linear and non-linear regression) suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. diff --git a/doc/pub/Regression/html/Regression-reveal.html b/doc/pub/Regression/html/Regression-reveal.html index c43d7c5b4..d5aaa438c 100644 --- a/doc/pub/Regression/html/Regression-reveal.html +++ b/doc/pub/Regression/html/Regression-reveal.html @@ -646,11 +646,11 @@ $$
-For a linear fit we don't need to invert a matrix!! - +For a linear fit we don't need to invert a matrix!! +Defining
$$
-\gamma = \sum_{i=0}^{n-1}\frac{n-1}{\sigma_i^2},
+\gamma = \sum_{i=0}^{n-1}\frac{1}{\sigma_i^2},
$$
@@ -678,7 +678,7 @@ $$
$$
-and show that
+we obtain
$$
\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
@@ -692,7 +692,7 @@ $$
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. +This approach (different linear and non-linear regression) suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below.
For a linear fit we don't need to invert a matrix!! +Defining $$ -\gamma = \sum_{i=0}^{n-1}\frac{n-1}{\sigma_i^2}, +\gamma = \sum_{i=0}^{n-1}\frac{1}{\sigma_i^2}, $$ @@ -637,7 +638,7 @@ $$ \gamma_{xy} = \sum_{i=0}^{n-1}\frac{y_ix_{i}}{\sigma_i^2}, $$ -and show that +we obtain $$ \beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, $$ @@ -647,7 +648,7 @@ $$ $$
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. +This approach (different linear and non-linear regression) suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. diff --git a/doc/pub/Regression/html/Regression.html b/doc/pub/Regression/html/Regression.html index 854e388de..7e1a6bbe5 100644 --- a/doc/pub/Regression/html/Regression.html +++ b/doc/pub/Regression/html/Regression.html @@ -621,8 +621,9 @@ $$
For a linear fit we don't need to invert a matrix!! +Defining $$ -\gamma = \sum_{i=0}^{n-1}\frac{n-1}{\sigma_i^2}, +\gamma = \sum_{i=0}^{n-1}\frac{1}{\sigma_i^2}, $$ @@ -642,7 +643,7 @@ $$ \gamma_{xy} = \sum_{i=0}^{n-1}\frac{y_ix_{i}}{\sigma_i^2}, $$ -and show that +we obtain $$ \beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, $$ @@ -652,7 +653,7 @@ $$ $$
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. +This approach (different linear and non-linear regression) suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below. diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb index e58efcf70..468d5907d 100644 --- a/doc/pub/Regression/ipynb/Regression.ipynb +++ b/doc/pub/Regression/ipynb/Regression.ipynb @@ -745,7 +745,8 @@ "source": [ "## The $\\chi^2$ function\n", "\n", - "For a linear fit we don't need to invert a matrix!!" + "For a linear fit we don't need to invert a matrix!! \n", + "Defining" ] }, { @@ -753,7 +754,7 @@ "metadata": {}, "source": [ "$$\n", - "\\gamma = \\sum_{i=0}^{n-1}\\frac{n-1}{\\sigma_i^2},\n", + "\\gamma = \\sum_{i=0}^{n-1}\\frac{1}{\\sigma_i^2},\n", "$$" ] }, @@ -842,7 +843,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "and show that" + "we obtain" ] }, { @@ -882,7 +883,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below.\n", + "This approach (different linear and non-linear regression) suffers often from both being underdetermined and overdetermined in the unknown coefficients $\\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed below. Or using Lasso and Ridge regression. See below.\n", "\n", "\n", "\n", diff --git a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz index 193056ca4..511f8a4c0 100644 Binary files a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz and b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz differ diff --git a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf index 17fab86a4..9a78bb544 100644 Binary files a/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf and b/doc/pub/Regression/pdf/Regression-beamer-handouts2x3.pdf differ diff --git a/doc/pub/Regression/pdf/Regression-beamer.pdf b/doc/pub/Regression/pdf/Regression-beamer.pdf index 1d9c4a61d..5d6b26295 100644 Binary files a/doc/pub/Regression/pdf/Regression-beamer.pdf and b/doc/pub/Regression/pdf/Regression-beamer.pdf differ diff --git a/doc/pub/Regression/pdf/Regression-minted.pdf b/doc/pub/Regression/pdf/Regression-minted.pdf index 2bfce68d9..d21efc3c3 100644 Binary files a/doc/pub/Regression/pdf/Regression-minted.pdf and b/doc/pub/Regression/pdf/Regression-minted.pdf differ diff --git a/doc/pub/Statistics/html/Statistics-bs.html b/doc/pub/Statistics/html/Statistics-bs.html index 5703370f9..a3574cb28 100644 --- a/doc/pub/Statistics/html/Statistics-bs.html +++ b/doc/pub/Statistics/html/Statistics-bs.html @@ -270,7 +270,7 @@ end of tocinfo --> -
- - - - -
- - -
-Another important distribution in science is the exponential distribution -$$ -\begin{equation*} -p(x) = \alpha\exp{-(\alpha x)}. -\end{equation*} -$$ -
- - -
-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 - -$$ -\begin{equation} -\langle h \rangle_X \equiv \int\! h(x)p(x)\,dx -\label{eq:expectation_value_of_h_wrt_p} -\end{equation} -$$ - -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 -$$ -\begin{equation*} -\langle x^n \rangle \equiv \int\! x^n p(x)\,dx -\end{equation*} -$$ -
- - -
-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 \) -$$ -\begin{equation*} -\langle x\rangle = \mu \equiv \int x p(x)dx, -\end{equation*} -$$ - -for a continuous distribution and -$$ -\begin{equation*} -\langle x\rangle = \mu \equiv \frac{1}{N}\sum_{i=1}^N x_i p(x_i), -\end{equation*} -$$ - -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) \). -
- - -
- -
-A special version of the moments is the set of central moments, the n-th central moment defined as -$$ -\begin{equation*} -\langle (x-\langle x\rangle )^n\rangle \equiv \int\! (x-\langle x\rangle)^n p(x)\,dx -\end{equation*} -$$ - -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) \) -$$ -\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*} -$$ - -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. -
- - -
- -
-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 | |
| 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 \) |
-
- - -
-With a PDF we can compute expectation values of selected quantities such as - -$$ -\begin{equation*} - \langle x^k\rangle=\frac{1}{N}\sum_{i=1}^{N}x_i^kp(x_i), -\end{equation*} -$$ - -if we have a discrete PDF or - -$$ -\begin{equation*} - \langle x^k\rangle=\int_a^b x^kp(x)dx, -\end{equation*} -$$ - -in the case of a continuous PDF. We have already defined the mean value \( \mu \) -and the variance \( \sigma^2 \). -
- - -
- -
-There are at least three PDFs which one may encounter. These are the - -
-Uniform distribution -$$ -\begin{equation*} -p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x), -\end{equation*} -$$ - -yielding probabilities different from zero in the interval \( [a,b] \). - -
-The exponential distribution -$$ -\begin{equation*} -p(x)=\alpha \exp{(-\alpha x)}, -\end{equation*} -$$ - -yielding probabilities different from zero in the interval \( [0,\infty) \) and with mean value -$$ -\begin{equation*} -\mu = \int_0^{\infty}xp(x)dx=\int_0^{\infty}x\alpha \exp{(-\alpha x)}dx=\frac{1}{\alpha}, -\end{equation*} -$$ -
- - -
-Finally, we have the so-called univariate normal distribution, or just the normal distribution -$$ -\begin{equation*} -p(x)=\frac{1}{b\sqrt{2\pi}}\exp{\left(-\frac{(x-a)^2}{2b^2}\right)} -\end{equation*} -$$ - -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 -$$ -\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*} -$$ - -which becomes with a suitable change of variables -$$ -\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*} -$$ -
- - -
-Similarly, the variance becomes -$$ -\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*} -$$ - -and inserting the mean value and performing a variable change we obtain - -$$ -\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*} -$$ - -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 - -$$ -\begin{equation} - p(x)=\frac{1}{\sqrt{2\pi}}\exp{\left(-\frac{x^2}{2}\right)}. -\label{_auto1} -\end{equation} -$$ -
- - -
- -
-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 - -$$ -\begin{equation*} -P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x\exp{\left(-\frac{t^2}{2}\right)}dt, -\end{equation*} -$$ - -which is difficult to evaluate in a quick way. -
- - -
- -
-Some other PDFs which one encounters often in the natural sciences are the binomial distribution -$$ -\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*} -$$ - -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 - -
- - -
- -
-In order to compute the mean and variance we need to recall Newton's binomial -formula -$$ -\begin{equation*} - (a+b)^m=\sum_{n=0}^m \left(\begin{array}{c} m \\ n\end{array}\right)a^nb^{m-n}, -\end{equation*} -$$ - -which can be used to show that - -$$ -\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*} -$$ - -the PDF is normalized to one. -The mean value is -$$ -\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*} -$$ - -resulting in -$$ -\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*} -$$ - -which we rewrite as - -$$ -\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*} -$$ -
- - -
- -
-Another important distribution with discrete stochastic variables \( x \) is -the Poisson model, which resembles the exponential distribution and reads -$$ -\begin{equation*} - p(x) = \frac{\lambda^x}{x!} e^{-\lambda} \hspace{0.5cm}x=0,1,\dots,;\lambda > 0. -\end{equation*} -$$ - -In this case both the mean value and the variance are easier to calculate, - -$$ -\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*} -$$ - -and the variance is \( \sigma^2=\lambda \). -
- - -
-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 - -$$ -\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*} -$$ -
- - -
-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 - -$$ -\begin{align} -\mathrm{Cov}(X_i,\,X_j) & = \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle -\label{_auto2}\\ -&=\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} -$$ - -with -$$ -\begin{equation*} -\langle x_i\rangle = -\int\cdots\int x_i P(x_1,\dots,x_n)\,dx_1\dots dx_n. -\end{equation*} -$$ -
- - -
-If we consider the above covariance as a matrix -$$ -C_{ij} =\mathrm{Cov}(X_i,\,X_j), -$$ - -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. -
- - -
- - -
# 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)
-- - -
- -
-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 \)) -$$ -\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*} -$$ - -If \( X_i \) and \( X_j \) are independent, we get -$$ -\langle x_i x_j\rangle = -\langle x_i\rangle\langle x_j\rangle=\mathrm{Cov}(X_i, X_j) = 0\ \ (i\neq j). -$$ -
- - -
- -
-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 -$$ -\begin{equation*} -\{x_1, x_2,\dots\,x_k,\dots\}. -\end{equation*} -$$ -
- - -
-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} \). -
- - -
-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 -$$ -\begin{equation*} -\langle x_{\alpha} \rangle \equiv \frac{1}{n}\sum_{k=1}^n x_{\alpha,k}. -\end{equation*} -$$ - -The sample variance is: -$$ -\begin{equation*} -\mathrm{Var}(x) \equiv \frac{1}{n}\sum_{k=1}^n (x_{\alpha,k} - \langle x_{\alpha} \rangle)^2, -\end{equation*} -$$ - -with its square root being the standard deviation of the sample. -
- - -
-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 -$$ -\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} -$$ - -where the last sums end at \( m \) and \( n \). -The total variance is -$$ -\begin{equation*} -\sigma^2_m= \frac{1}{mn^2}\sum_{\alpha=1}^m(\langle x_{\alpha} \rangle-\langle X_m \rangle)^2, -\end{equation*} -$$ - -which we rewrite as -$$ -\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} -$$ -
- - -
- -
-We define also the sample variance \( \sigma^2 \) of all \( mn \) individual experiments as -$$ -\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} -$$ - -
-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) \). -
- - -
- -
-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 -$$ -\begin{equation*} - \sigma_m= -\frac{\sigma}{\sqrt{m}}. -\end{equation*} -$$ - -
-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. -
- - -
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \) - -$$ -\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*} -$$ - -
-We can then use Eq. \eqref{eq:exptvariance} -$$ -\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*} -$$ - -and rewrite it as -$$ -\begin{equation*} -\sigma^2_m=\frac{\sigma^2}{n}+\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k < l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle), -\end{equation*} -$$ - -where the first term is the sample variance of all \( mn \) experiments divided by \( n \) -and the last term is nothing but the covariance which arises when \( k\ne l \). -
- - -
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \) - -
-If the -observables are uncorrelated, then the covariance is zero and we obtain a total variance -which agrees with the central limit theorem. Correlations may often be present in our data set, resulting in a non-zero covariance. The first term is normally called the uncorrelated -contribution. -Computationally the uncorrelated first term is much easier to treat -efficiently than the second. -We just accumulate separately the values \( x^2 \) and \( x \) for every -measurement \( x \) we receive. The correlation term, though, has to be -calculated at the end of the experiment since we need all the -measurements to calculate the cross terms. Therefore, all measurements -have to be stored throughout the experiment. -
- - -
- -
-Let us analyze the problem by splitting up the correlation term into -partial sums of the form - -$$ -\begin{equation*} -f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle), -\end{equation*} -$$ - -The correlation term of the total variance can now be rewritten in terms of -\( f_d \) - -$$ -\begin{equation*} -\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k < l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle)= -\frac{2}{n}\sum_{d=1}^{n-1} f_d -\end{equation*} -$$ -
- - -
-The value of \( f_d \) reflects the correlation between measurements -separated by the distance \( d \) in the samples. Notice that for -\( d=0 \), \( f \) is just the sample variance, \( \sigma^2 \). If we divide \( f_d \) -by \( \sigma^2 \), we arrive at the so called autocorrelation function - -$$ -\begin{equation} -\kappa_d = \frac{f_d}{\sigma^2} -\label{eq:autocorrelformal} -\end{equation} -$$ - -which gives us a useful measure of the correlation pair correlation -starting always at \( 1 \) for \( d=0 \). -
- - -
- -
-The sample variance of the \( mn \) experiments can now be -written in terms of the autocorrelation function - -$$ -\begin{equation} -\sigma_m^2=\frac{\sigma^2}{n}+\frac{2}{n}\cdot\sigma^2\sum_{d=1}^{n-1} -\frac{f_d}{\sigma^2}=\left(1+2\sum_{d=1}^{n-1}\kappa_d\right)\frac{1}{n}\sigma^2=\frac{\tau}{n}\cdot\sigma^2 -\label{eq:error_estimate_corr_time} -\end{equation} -$$ - -and we see that \( \sigma_m \) can be expressed in terms of the -uncorrelated sample variance times a correction factor \( \tau \) which -accounts for the correlation between measurements. We call this -correction factor the autocorrelation time - -$$ -\begin{equation} -\tau = 1+2\sum_{d=1}^{n-1}\kappa_d -\label{eq:autocorrelation_time} -\end{equation} -$$ - - - -For a correlation free experiment, \( \tau \) -equals 1. -
- - -
-From the point of view of -Eq. \eqref{eq:error_estimate_corr_time} we can interpret a sequential -correlation as an effective reduction of the number of measurements by -a factor \( \tau \). The effective number of measurements becomes -$$ -\begin{equation*} -n_\mathrm{eff} = \frac{n}{\tau} -\end{equation*} -$$ - -To neglect the autocorrelation time \( \tau \) will always cause our -simple uncorrelated estimate of \( \sigma_m^2\approx \sigma^2/n \) to -be less than the true sample error. The estimate of the error will be -too "good". On the other hand, the calculation of the full -autocorrelation time poses an efficiency problem if the set of -measurements is very large. The solution to this problem is given by -more practically oriented methods like the blocking technique. - -
- - -
- - -
# Importing various packages
-from math import exp, sqrt
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Sample covariance, note the factor 1/(n-1)
-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-1.)
-
-n = 100
-x = np.random.normal(size=n)
-print(np.mean(x))
-y = 4+3*x+np.random.normal(size=n)
-print(np.mean(y))
-z = x**3+np.random.normal(size=n)
-print(np.mean(z))
-covxx = covariance(x,x,n)
-covyy = covariance(y,y,n)
-covzz = covariance(z,z,n)
-covxy = covariance(x,y,n)
-covxz = covariance(x,z,n)
-covyz = covariance(y,z,n)
-print(covxx,covyy, covzz)
-print(covxy,covxz, covyz)
-w = np.vstack((x, y, z))
-#print(w)
-c = np.cov(w)
-print(c)
-#eigen = np.zeros(n)
-Eigvals, Eigvecs = np.linalg.eig(c)
-print(Eigvals)
-- - -
- -
-Uniform deviates are just random numbers that lie within a specified range -(typically 0 to 1), with any one number in the range just as likely as any other. They -are, in other words, what you probably think random numbers are. However, -we want to distinguish uniform deviates from other sorts of random numbers, for -example numbers drawn from a normal (Gaussian) distribution of specified mean -and standard deviation. These other sorts of deviates are almost always generated by -performing appropriate operations on one or more uniform deviates, as we will see -in subsequent sections. So, a reliable source of random uniform deviates, the subject -of this section, is an essential building block for any sort of stochastic modeling -or Monte Carlo computer work. -
- - -
- -
-A disclaimer is however appropriate. It should be fairly obvious that -something as deterministic as a computer cannot generate purely random numbers. - -
-Numbers generated by any of the standard algorithms are in reality pseudo random -numbers, hopefully abiding to the following criteria: - -
- - -
- The most common random number generators are based on so-called -Linear congruential relations of the type - -$$ -\begin{equation*} - N_i=(aN_{i-1}+c) \mathrm{MOD} (M), -\end{equation*} -$$ - -which yield a number in the interval [0,1] through - -$$ -\begin{equation*} - x_i=N_i/M -\end{equation*} -$$ - -
-The number -\( M \) is called the period and it should be as large as possible - and -\( N_0 \) is the starting value, or seed. The function \( \mathrm{MOD} \) means the remainder, -that is if we were to evaluate \( (13)\mathrm{MOD}(9) \), the outcome is the remainder -of the division \( 13/9 \), namely \( 4 \). -
- - -
- -
-The problem with such generators is that their outputs are periodic; -they -will start to repeat themselves with a period that is at most \( M \). If however -the parameters \( a \) and \( c \) are badly chosen, the period may be even shorter. - -
-Consider the following example - -$$ -\begin{equation*} - N_i=(6N_{i-1}+7) \mathrm{MOD} (5), -\end{equation*} -$$ - -with a seed \( N_0=2 \). This generator produces the sequence -\( 4,1,3,0,2,4,1,3,0,2,...\dots \), i.e., a sequence with period \( 5 \). -However, increasing \( M \) may not guarantee a larger period as the following -example shows - -$$ -\begin{equation*} - N_i=(27N_{i-1}+11) \mathrm{MOD} (54), -\end{equation*} -$$ - -which still, with \( N_0=2 \), results in \( 11,38,11,38,11,38,\dots \), a period of -just \( 2 \). -
- - -
-Typical periods for the random generators provided in the program library -are of the order of \( \sim 10^9 \) or larger. Other random number generators which have -become increasingly popular are so-called shift-register generators. -In these generators each successive number depends on many preceding -values (rather than the last values as in the linear congruential -generator). -For example, you could make a shift register generator whose $l$th -number is the sum of the $l-i$th and $l-j$th values with modulo \( M \), -$$ -\begin{equation*} - N_l=(aN_{l-i}+cN_{l-j})\mathrm{MOD}(M). -\end{equation*} -$$ -
- - -
-Such a generator again produces a sequence of pseudorandom numbers -but this time with a period much larger than \( M \). -It is also possible to construct more elaborate algorithms by including -more than two past terms in the sum of each iteration. -One example is the generator of Marsaglia and Zaman -which consists of two congruential relations - -$$ -\begin{equation} - N_l=(N_{l-3}-N_{l-1})\mathrm{MOD}(2^{31}-69), -\label{eq:mz1} -\end{equation} -$$ - -followed by -$$ -\begin{equation} - N_l=(69069N_{l-1}+1013904243)\mathrm{MOD}(2^{32}), -\label{eq:mz2} -\end{equation} -$$ - -which according to the authors has a period larger than \( 2^{94} \). -
- - -
-Instead of using modular addition, we could use the bitwise -exclusive-OR (\( \oplus \)) operation so that - -$$ -\begin{equation*} - N_l=(N_{l-i})\oplus (N_{l-j}) -\end{equation*} -$$ - -where the bitwise action of \( \oplus \) means that if \( N_{l-i}=N_{l-j} \) the result is -\( 0 \) whereas if \( N_{l-i}\ne N_{l-j} \) the result is -\( 1 \). As an example, consider the case where \( N_{l-i}=6 \) and \( N_{l-j}=11 \). The first -one has a bit representation (using 4 bits only) which reads \( 0110 \) whereas the -second number is \( 1011 \). Employing the \( \oplus \) operator yields -\( 1101 \), or \( 2^3+2^2+2^0=13 \). - -
-In Fortran90, the bitwise \( \oplus \) operation is coded through the intrinsic -function \( \mathrm{IEOR}(m,n) \) where \( m \) and \( n \) are the input numbers, while in \( C \) -it is given by \( m\wedge n \). -
- - -
- -
-We show here how the linear congruential algorithm can be implemented, namely -$$ -\begin{equation*} - N_i=(aN_{i-1}) \mathrm{MOD} (M). -\end{equation*} -$$ - -However, since \( a \) and \( N_{i-1} \) are integers and their multiplication -could become greater than the standard 32 bit integer, there is a trick via -Schrage's algorithm which approximates the multiplication -of large integers through the factorization -$$ -\begin{equation*} - M=aq+r, -\end{equation*} -$$ - -where we have defined - -$$ -\begin{equation*} - q=[M/a], -\end{equation*} -$$ - -and -$$ -\begin{equation*} - r = M\hspace{0.1cm}\mathrm{MOD} \hspace{0.1cm}a. -\end{equation*} -$$ - -where the brackets denote integer division. In the code below the numbers -\( q \) and \( r \) are chosen so that \( r < q \). -
- - -
- -
-To see how this works we note first that -$$ -\begin{equation} -(aN_{i-1}) \mathrm{MOD} (M)= (aN_{i-1}-[N_{i-1}/q]M)\mathrm{MOD} (M), -\label{eq:rntrick1} -\end{equation} -$$ - -since we can add or subtract any integer multiple of \( M \) from \( aN_{i-1} \). -The last term \( [N_{i-1}/q]M\mathrm{MOD}(M) \) is zero since the integer division -\( [N_{i-1}/q] \) just yields a constant which is multiplied with \( M \). -
- - -
-We can now rewrite Eq. \eqref{eq:rntrick1} as - -$$ -\begin{equation} -(aN_{i-1}) \mathrm{MOD} (M)= (aN_{i-1}-[N_{i-1}/q](aq+r))\mathrm{MOD} (M), -\label{eq:rntrick2} -\end{equation} -$$ - -which results -in - -$$ -\begin{equation} -(aN_{i-1}) \mathrm{MOD} (M)= \left(a(N_{i-1}-[N_{i-1}/q]q)-[N_{i-1}/q]r)\right)\mathrm{MOD} (M), -\label{eq:rntrick3} -\end{equation} -$$ - -yielding -$$ -\begin{equation} -(aN_{i-1}) \mathrm{MOD} (M)= \left(a(N_{i-1}\mathrm{MOD} (q)) -[N_{i-1}/q]r)\right)\mathrm{MOD} (M). -\label{eq:rntrick4} -\end{equation} -$$ -
- - -
-The term \( [N_{i-1}/q]r \) is always smaller or equal \( N_{i-1}(r/q) \) and with \( r < q \) we obtain always a -number smaller than \( N_{i-1} \), which is smaller than \( M \). -And since the number \( N_{i-1}\mathrm{MOD} (q) \) is between zero and \( q-1 \) then -\( a(N_{i-1}\mathrm{MOD} (q)) < aq \). Combined with our definition of \( q=[M/a] \) ensures that -this term is also smaller than \( M \) meaning that both terms fit into a -32-bit signed integer. None of these two terms can be negative, but their difference could. -The algorithm below adds \( M \) if their difference is negative. -Note that the program uses the bitwise \( \oplus \) operator to generate -the starting point for each generation of a random number. The period -of \( ran0 \) is \( \sim 2.1\times 10^{9} \). A special feature of this -algorithm is that is should never be called with the initial seed -set to \( 0 \). -
- - -
- -
- - -
/*
- ** The function
- ** ran0()
- ** is an "Minimal" random number generator of Park and Miller
- ** Set or reset the input value
- ** idum to any integer value (except the unlikely value MASK)
- ** to initialize the sequence; idum must not be altered between
- ** calls for sucessive deviates in a sequence.
- ** The function returns a uniform deviate between 0.0 and 1.0.
- */
-double ran0(long &idum)
-{
- const int a = 16807, m = 2147483647, q = 127773;
- const int r = 2836, MASK = 123459876;
- const double am = 1./m;
- long k;
- double ans;
- idum ^= MASK;
- k = (*idum)/q;
- idum = a*(idum - k*q) - r*k;
- // add m if negative difference
- if(idum < 0) idum += m;
- ans=am*(idum);
- idum ^= MASK;
- return ans;
-} // End: function ran0()
--
- - -
- -
-As mentioned previously, the underlying PDF for the generation of -random numbers is the uniform distribution, meaning that the -probability for finding a number \( x \) in the interval [0,1] is \( p(x)=1 \). - -
-A random number generator should produce numbers which are uniformly distributed -in this interval. The table shows the distribution of \( N=10000 \) random -numbers generated by the functions in the program library. -We note in this table that the number of points in the various -intervals \( 0.0-0.1 \), \( 0.1-0.2 \) etc are fairly close to \( 1000 \), with some minor -deviations. - -
-Two additional measures are the standard deviation \( \sigma \) and the mean -\( \mu=\langle x\rangle \). -
- - -
-For the uniform distribution, the mean value \( \mu \) is then - -$$ -\begin{equation*} - \mu=\langle x\rangle=\frac{1}{2} -\end{equation*} -$$ - -while the standard deviation is - -$$ -\begin{equation*} - \sigma=\sqrt{\langle x^2\rangle-\mu^2}=\frac{1}{\sqrt{12}}=0.2886. -\end{equation*} -$$ -
- - -
-The various random number generators produce results which agree rather well with -these limiting values. - -
- -
| \( x \)-bin | ran0 | ran1 | ran2 | ran3 |
| 0.0-0.1 | 1013 | 991 | 938 | 1047 |
| 0.1-0.2 | 1002 | 1009 | 1040 | 1030 |
| 0.2-0.3 | 989 | 999 | 1030 | 993 |
| 0.3-0.4 | 939 | 960 | 1023 | 937 |
| 0.4-0.5 | 1038 | 1001 | 1002 | 992 |
| 0.5-0.6 | 1037 | 1047 | 1009 | 1009 |
| 0.6-0.7 | 1005 | 989 | 1003 | 989 |
| 0.7-0.8 | 986 | 962 | 985 | 954 |
| 0.8-0.9 | 1000 | 1027 | 1009 | 1023 |
| 0.9-1.0 | 991 | 1015 | 961 | 1026 |
| \( \mu \) | 0.4997 | 0.5018 | 0.4992 | 0.4990 |
| \( \sigma \) | 0.2882 | 0.2892 | 0.2861 | 0.2915 |
-
- - -
-The following simple Python code plots the distribution of the produced random numbers using the linear congruential RNG employed by Python. The trend displayed in the previous table is seen rather clearly. -
- - -
--
- - -
-Since our random numbers, which are typically generated via a linear congruential algorithm, -are never fully independent, we can then define -an important test which measures the degree of correlation, namely the so-called -auto-correlation function defined previously, see again Eq. \eqref{eq:autocorrelformal}. -We rewrite it here as -$$ -\begin{equation*} - C_k=\frac{f_d} - {\sigma^2}, -\end{equation*} -$$ - -with \( C_0=1 \). Recall that -\( \sigma^2=\langle x_i^2\rangle-\langle x_i\rangle^2 \) and that -$$ -\begin{equation*} -f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle), -\end{equation*} -$$ - -
-The non-vanishing of \( C_k \) for \( k\ne 0 \) means that the random -numbers are not independent. The independence of the random numbers is crucial -in the evaluation of other expectation values. If they are not independent, our -assumption for approximating \( \sigma_N \) is no longer valid. - -
-
- - -
- - -
# Importing various packages
-from math import exp, sqrt
-from random import random, seed
-import numpy as np
-import matplotlib.pyplot as plt
-
-def autocovariance(x, n, k, mean_x):
- sum = 0.0
- for i in range(0, n-k):
- sum += (x[(i+k)]-mean_x)*(x[i]-mean_x)
- return sum/n
-
-n = 1000
-x=np.random.normal(size=n)
-autocor = np.zeros(n)
-figaxis = np.zeros(n)
-mean_x=np.mean(x)
-var_x = np.var(x)
-print(mean_x, var_x)
-for i in range (0, n):
- figaxis[i] = i
- autocor[i]=(autocovariance(x, n, i, mean_x))/var_x
-
-plt.plot(figaxis, autocor, "r-")
-plt.axis([0,n,-0.1, 1.0])
-plt.xlabel(r'$i$')
-plt.ylabel(r'$\gamma_i$')
-plt.title(r'Autocorrelation function')
-plt.show()
--As can be seen from the plot, the first point gives back the variance and a value of one. -For the remaining values we notice that there are still non-zero values for the auto-correlation function. - -
- - -
-The program here computes the correlation function for one of the standard functions included with the c++ compiler. -
- - -
// This function computes the autocorrelation function for
-// the standard c++ random number generator
-
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <cmath>
-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
--
- - -
- -
- - -
-The following part of a c++ code (from project 4) sets up the uniform distribution for \( x\in [0,1] \). -
- - -
/*
-
-// You need this
-#include <random>
-
-// 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<double> RandomNumberGenerator(0.0,1.0);
-
-// Now use the RNG
-int ix = (int) (RandomNumberGenerator(gen)*NSpins);
--
- - -
- -
-
- - -
- -
- - -
- - -
// This function computes the autocorrelation function for
-// the Mersenne random number generator with a uniform distribution
-#include <iostream>
-#include <fstream>
-#include <iomanip>
-#include <cstdlib>
-#include <random>
-#include <armadillo>
-#include <string>
-#include <cmath>
-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<vec>(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<double> 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<vec>(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
-- - -
- -
- - -
- -
- - -
- -
- - -
- -
- - -
- - -
# 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)
-- - -
-The Bootstrap resampling method is also very popular. It is very simple: - -
-This method is very useful for small ensembles of data points. - -
- - -
- - -
-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}} \) . - -
- - -
- - -
# 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)
-- - -
- - -
# 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)
-- diff --git a/doc/pub/Statistics/html/Statistics-reveal.html b/doc/pub/Statistics/html/Statistics-reveal.html index b5b596ff1..c80275c6a 100644 --- a/doc/pub/Statistics/html/Statistics-reveal.html +++ b/doc/pub/Statistics/html/Statistics-reveal.html @@ -163,7 +163,7 @@ MathJax.Hub.Config({
-
@@ -176,6 +176,7 @@ MathJax.Hub.Config({
Things to add
Add general statistic elements (probability theory mainly), assumed knowledge
+Repeat basic definitions.
| Discrete PDF | Continuous PDF | |
|---|---|---|
| Discrete PDF | Continuous PDF | |
| 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 \) |
| 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 \) |
-
@@ -741,16 +742,16 @@ while \( P(x) \) is the cumulative probability.
| Discrete PDF | Continuous PDF | |
|---|---|---|
| Discrete PDF | Continuous PDF | |
| 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 \) |
| 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 \) |
-
@@ -746,16 +747,16 @@ while \( P(x) \) is the cumulative probability.
| Discrete PDF | Continuous PDF | |
|---|---|---|
| Discrete PDF | Continuous PDF | |
| 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 \) |
| 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 \) |