FYS-STK4155/doc/LectureNotes/_build/jupyter_execute/statistics.ipynb

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

\{2,3,4,5,6,7,8,9,10,11,12\}.

These values are called the domain. To this domain we have the corresponding probabilities

\{1/36,2/36/,3/36,4/36,5/36,6/36,5/36,4/36,3/36,2/36,1/36\}.

Tossing the dice

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

\{2,3,4,5,6,7,8,9,10,11,12\},

appear in a random order. After 11 throws the results may look like

\{10,8,6,3,6,9,11,8,12,4,5\}.

Stochastic variables

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

p(x) = \mathrm{Prob}(X=x).

Stochastic variables and the main concepts, the continuous case

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

\mathrm{Prob}(a\leq X\leq b) = \int_a^b p(x)dx.

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.

The cumulative probability

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

P(x)=\mathrm{Prob(}X\leq x\mathrm{)} =
\int_{-\infty}^x p(x^{\prime})dx^{\prime}.

The relation between a CDF and its corresponding PDF is then

p(x) = \frac{d}{dx}P(x).

Properties of PDFs

There are two properties that all PDFs must satisfy. The first one is positivity (assuming that the PDF is normalized)

0 \leq p(x) \leq 1.

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

\begin{align*}
\sum_{x_i\in\mathbb D} p(x_i) & =  1,\\
\int_{x\in\mathbb D} p(x)\,dx & =  1.
\end{align*}

Important distributions, the uniform distribution

The first one is the most basic PDF; namely the uniform distribution

\begin{equation}
p(x) = \frac{1}{b-a}\theta(x-a)\theta(b-x).
\label{eq:unifromPDF} \tag{1}
\end{equation}

For a=0 and b=1 we have

\begin{array}{ll}
p(x)dx = dx & \in [0,1].
\end{array}

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.

Gaussian distribution

The second one is the Gaussian Distribution

p(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp{(-\frac{(x-\mu)^2}{2\sigma^2})},

with mean value \mu and standard deviation \sigma. If \mu=0 and \sigma=1, it is normally called the standard normal distribution

p(x) = \frac{1}{\sqrt{2\pi}} \exp{(-\frac{x^2}{2})},

The following simple Python code plots the above distribution for different values of \mu and \sigma.

In [1]:

%matplotlib inline

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()

Exponential distribution

Another important distribution in science is the exponential distribution

p(x) = \alpha\exp{-(\alpha x)}.

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

\begin{equation}
\langle h \rangle_X \equiv \int\! h(x)p(x)\,dx
\label{eq:expectation_value_of_h_wrt_p} \tag{2}
\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

\langle x^n \rangle \equiv \int\! x^n p(x)\,dx

Stochastic variables and the main concepts, mean values

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

\langle x\rangle  = \mu \equiv \int x p(x)dx,

for a continuous distribution and

\langle x\rangle  = \mu \equiv \sum_{i=1}^N x_i p(x_i),

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).

Stochastic variables and the main concepts, central moments, the variance

A special version of the moments is the set of central moments, the n-th central moment defined as

\langle (x-\langle x\rangle )^n\rangle  \equiv \int\! (x-\langle x\rangle)^n p(x)\,dx

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 - 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.

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
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$

Probability Distribution Functions

With a PDF we can compute expectation values of selected quantities such as

\langle x^k\rangle=\sum_{i=1}^{N}x_i^kp(x_i),

if we have a discrete PDF or

\langle x^k\rangle=\int_a^b x^kp(x)dx,

in the case of a continuous PDF. We have already defined the mean value \mu and the variance \sigma^2.

The three famous Probability Distribution Functions

There are at least three PDFs which one may encounter. These are the

Uniform distribution

p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x),

yielding probabilities different from zero in the interval [a,b].

The exponential distribution

p(x)=\alpha \exp{(-\alpha x)},

yielding probabilities different from zero in the interval [0,\infty) and with mean value

\mu = \int_0^{\infty}xp(x)dx=\int_0^{\infty}x\alpha \exp{(-\alpha x)}dx=\frac{1}{\alpha},

with variance

\sigma^2=\int_0^{\infty}x^2p(x)dx-\mu^2 = \frac{1}{\alpha^2}.

Probability Distribution Functions, the normal distribution

Finally, we have the so-called univariate normal distribution, or just the normal distribution

p(x)=\frac{1}{b\sqrt{2\pi}}\exp{\left(-\frac{(x-a)^2}{2b^2}\right)}

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

\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,

which becomes with a suitable change of variables

\mu =\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(a+b\sqrt{2}y)\exp{-y^2}dy=a.

Probability Distribution Functions, the normal distribution

Similarly, the variance becomes

\sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}(x-\mu)^2 \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx,

and inserting the mean value and performing a variable change we obtain

\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,

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} \tag{3}
\end{equation}

Probability Distribution Functions, the cumulative distribution

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

P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x\exp{\left(-\frac{t^2}{2}\right)}dt,

which is difficult to evaluate in a quick way.

Probability Distribution Functions, other important distribution

Some other PDFs which one encounters often in the natural sciences are the binomial distribution

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,

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.

Probability Distribution Functions, the binomial distribution

In order to compute the mean and variance we need to recall Newton's binomial formula

(a+b)^m=\sum_{n=0}^m \left(\begin{array}{c} m \\ n\end{array}\right)a^nb^{m-n},

which can be used to show that

\sum_{x=0}^n\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} = (y+1-y)^n = 1,

the PDF is normalized to one. The mean value is

\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},

resulting in

\mu = 
\sum_{x=0}^n x\frac{(n-1)!}{(x-1)!(n-1-(x-1))!}y^{x-1}(1-y)^{n-1-(x-1)},

which we rewrite as

\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.

The variance is slightly trickier to get. It reads \sigma^2=ny(1-y).

Probability Distribution Functions, Poisson's distribution

Another important distribution with discrete stochastic variables x is
the Poisson model, which resembles the exponential distribution and reads

p(x) = \frac{\lambda^x}{x!} e^{-\lambda} \hspace{0.5cm}x=0,1,\dots,;\lambda > 0.

In this case both the mean value and the variance are easier to calculate,

\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,

and the variance is \sigma^2=\lambda.

Probability Distribution Functions, Poisson's distribution

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

\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}.

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

\begin{equation}
\mathrm{Cov}(X_i,\,X_j)  = \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle 
\label{_auto2} \tag{4}
\end{equation}

\begin{equation} 
=\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} \tag{5}
\end{equation}

with

\langle x_i\rangle =
\int\cdots\int x_i P(x_1,\dots,x_n)\,dx_1\dots dx_n.

Meet the covariance in matrix disguise

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.

Covariance

In [2]:

# 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)

4.323291478597321
[[4.06126507e+00 5.22717936e+00 3.58228342e-01 1.09103481e+01
  4.42847770e+00 5.02161783e+00 8.06412177e-03 6.67407338e+00
  1.12366979e+01 6.04205220e+00]
 [5.22717936e+00 6.72780613e+00 4.61069091e-01 1.40425078e+01
  5.69981195e+00 6.46323168e+00 1.03791824e-02 8.59007674e+00
  1.44625466e+01 7.77661393e+00]
 [3.58228342e-01 4.61069091e-01 3.15979239e-02 9.62359224e-01
  3.90618734e-01 4.42937310e-01 7.11304709e-04 5.88693966e-01
  9.91145266e-01 5.32945844e-01]
 [1.09103481e+01 1.40425078e+01 9.62359224e-01 2.93100040e+01
  1.18968431e+01 1.34902789e+01 2.16637855e-02 1.79295029e+01
  3.01867234e+01 1.62316154e+01]
 [4.42847770e+00 5.69981195e+00 3.90618734e-01 1.18968431e+01
  4.82889306e+00 5.47566390e+00 8.79326583e-03 7.27753165e+00
  1.22527008e+01 6.58836420e+00]
 [5.02161783e+00 6.46323168e+00 4.42937310e-01 1.34902789e+01
  5.47566390e+00 6.20906175e+00 9.97101567e-03 8.25226753e+00
  1.38937995e+01 7.47079457e+00]
 [8.06412177e-03 1.03791824e-02 7.11304709e-04 2.16637855e-02
  8.79326583e-03 9.97101567e-03 1.60122668e-05 1.32521615e-02
  2.23117916e-02 1.19972087e-02]
 [6.67407338e+00 8.59007674e+00 5.88693966e-01 1.79295029e+01
  7.27753165e+00 8.25226753e+00 1.32521615e-02 1.09678277e+01
  1.84658093e+01 9.92919670e+00]
 [1.12366979e+01 1.44625466e+01 9.91145266e-01 3.01867234e+01
  1.22527008e+01 1.38937995e+01 2.23117916e-02 1.84658093e+01
  3.10896672e+01 1.67171347e+01]
 [6.04205220e+00 7.77661393e+00 5.32945844e-01 1.62316154e+01
  6.58836420e+00 7.47079457e+00 1.19972087e-02 9.92919670e+00
  1.67171347e+01 8.98892195e+00]]

Meet the covariance, uncorrelated events

Consider the stochastic variables X_i and X_j, (i\neq j). We have

\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*}

If X_i and X_j are independent (assuming i \neq j), we have that

Warning:

Output truncated. This notebook contains too many cells to display efficiently.