======= Elements of Probability Theory and Statistical Data Analysis ======= ===== Domains and probabilities ===== 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 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 _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)_. === Stochastic variables and the main concepts, the discrete case === 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 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. 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 === Properties of PDFs === 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 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. 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 Another important distribution in science is the exponential distribution !bt \begin{equation*} p(x) = \alpha\exp{-(\alpha x)}. \end{equation*} !et === Expectation values === 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 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)$. 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. === Probability Distribution Functions === 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$ | |--------------------------------------------------------------------------------------------------------------------------------------| 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$. 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 with variance !bt \begin{equation*} \sigma^2=\int_0^{\infty}x^2p(x)dx-\mu^2 = \frac{1}{\alpha^2}. \end{equation*} !et 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 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 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. 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$. 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 The variance is slightly trickier to get. It reads $\sigma^2=ny(1-y)$. 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$. 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 === Meet the covariance! === 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 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. !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 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 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 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 !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*. 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 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)$. === Numerical experiments and the covariance, central limit theorem === 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. 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