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

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

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

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.

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

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.

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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
      3 import numpy as np
      4 from math import acos, exp, sqrt

File ~/miniforge3/lib/python3.9/site-packages/IPython/core/interactiveshell.py:2432, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2430     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2431 with self.builtin_trap:
-> 2432     result = fn(*args, **kwargs)
   2434 # The code below prevents the output from being displayed
   2435 # when using magics with decorator @output_can_be_silenced
   2436 # when the last Python token in the expression is a ';'.
   2437 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/miniforge3/lib/python3.9/site-packages/IPython/core/magics/pylab.py:99, in PylabMagics.matplotlib(self, line)
     97     print("Available matplotlib backends: %s" % backends_list)
     98 else:
---> 99     gui, backend = self.shell.enable_matplotlib(args.gui.lower() if isinstance(args.gui, str) else args.gui)
    100     self._show_matplotlib_backend(args.gui, backend)

File ~/miniforge3/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3606, in InteractiveShell.enable_matplotlib(self, gui)
   3585 def enable_matplotlib(self, gui=None):
   3586     """Enable interactive matplotlib and inline figure support.
   3587 
   3588     This takes the following steps:
   (...)
   3604         display figures inline.
   3605     """
-> 3606     from matplotlib_inline.backend_inline import configure_inline_support
   3608     from IPython.core import pylabtools as pt
   3609     gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)

File ~/miniforge3/lib/python3.9/site-packages/matplotlib_inline/__init__.py:1
----> 1 from . import backend_inline, config  # noqa
      2 __version__ = "0.1.6"  # noqa

File ~/miniforge3/lib/python3.9/site-packages/matplotlib_inline/backend_inline.py:6
      1 """A matplotlib backend for publishing figures via display_data"""
      3 # Copyright (c) IPython Development Team.
      4 # Distributed under the terms of the BSD 3-Clause License.
----> 6 import matplotlib
      7 from matplotlib import colors
      8 from matplotlib.backends import backend_agg

ModuleNotFoundError: No module named 'matplotlib'

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

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

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$

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.

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

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.

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}

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.

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.

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

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.

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.

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.

Warning:

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