Statistics

Abstract

Statistics is a fundamental tool to study physical phenomena in any system, be it complex or not. However, it is particularly important in the study of complex systems, for reasons that will be explained in this chapter in detail.

Appendices

Appendix 1: The Fourier Transform

The Fourier representation has a very large number of applications in Mathematics, Physics and Engineering, since it permits to express an arbitrary function as a linear combination of periodic functions. Truncation of the Fourier representation, for example, is the basis of signal filtering. Many other manipulations (modulation, integration, etc.) form the basis of modern communications. In Physics, use of the Fourier representation is essential for the understanding of resonances, oscillations, wave propagation, transmission and absorption, fluid and plasma turbulence and many other processes. In Mathematics, Fourier transforms are often used to simplify the solution of differential equations, among many other problems.

The Fourier transform of a (maybe complex) function f(x) is defined as⁵⁸:

$$\displaystyle \begin{aligned} \mathrm{F}[\,f(x)] := \hat f(k) = \int_{-\infty}^\infty dx~ f(x) \exp(\i kx), \end{aligned} $$

(2.123)

where $\i = \sqrt {-1}$, k is known as the wave number.⁵⁹ Note that, if f(x) is real, then $\hat f(k)$ is Hermitian, meaning that $\hat f(-k) = \hat f^*(k)$, where the asterisk represents the complex conjugate.

The inverse of the Fourier transform is then obtained as,

$$\displaystyle \begin{aligned} f(x) = \mathrm{F}^{-1}[\hat f (k)]: = \frac{1}{2\pi}\int_{-\infty}^\infty dk~ \hat f(k) \exp(-\i kx ). \end{aligned} $$

(2.124)

This expression is also referred to as the Fourier representation of f(x).

Some common Fourier transform pairs are collected in Table 2.7. However, it must be noted that the existence of both the Fourier and inverse transforms is not always guaranteed. A sufficient condition (but not necessary) is that both f(x) and $\hat f(k)$ be (Lebesgue)-integrable. That is,

$$\displaystyle \begin{aligned} \int_{-\infty}^\infty dx~|\,f(x)| < \infty,~~~~\int_{-\infty}^\infty dk~|\,\hat f(k)| < \infty, \end{aligned} $$

(2.125)

which requires that | f(x)|→ 0 faster than |x|⁻¹ for x →±∞, and that $|\,\hat f(k)|\rightarrow 0$ faster than |k|⁻¹ for k →±∞. But Fourier transforms may exist for some functions that violate this condition.

Table 2.7 Some useful Fourier transforms

Fourier transforms have many interesting properties. In particular, they are linear operations. That is, the Fourier transform of the sum of any two functions is equal to the sum of their Fourier transforms.

Secondly, the Fourier transform of the n -th derivative of a function becomes particularly simple:

$$\displaystyle \begin{aligned} \mathrm{F}\left[\frac{d^n f}{dx^n}\right] = (-\i k)^n \hat f(k), \end{aligned} $$

(2.126)

as trivially follows from differentiating Eq. 2.124.

Another interesting property has to do with the translation of the argument, that translates into a multiplication by a phase in Fourier space:

$$\displaystyle \begin{aligned} \mathrm{F}[\,f(x\pm x_0)] = \exp(\mp\i kx_0)\hat f(k). \end{aligned} $$

(2.127)

Particularly useful properties are the following two theorems [23]. One is Parseval’s theorem,

$$\displaystyle \begin{aligned} \int_{-\infty}^\infty dx~ f(x) g^*(x) =\frac{1}{2\pi} \int_{-\infty}^\infty dk~ \hat f(k) \hat g^*(k) \end{aligned} $$

(2.128)

where the asterisk denotes complex conjugation. When f = g, this becomes,

$$\displaystyle \begin{aligned} \int_{-\infty}^\infty dx~ |\,f(x)|{}^2 =\frac{1}{2\pi} \int_{-\infty}^\infty dk~ |\,\hat f(k)|{}^2. \end{aligned} $$

(2.129)

$|\,\hat f(k)|{ }^2$ is usually known as the power spectrum of f(x).⁶⁰

The second important theorem is the convolution theorem. It states that the convolution of two functions, defined as:

$$\displaystyle \begin{aligned} c_{[\,f,g]}(x) = \int_{-\infty}^\infty f(x')g(x-x') dx' \end{aligned} $$

(2.130)

satisfies that its Fourier transform is given by:

$$\displaystyle \begin{aligned} \hat c_{[\,f,g]}(k) = \hat f(k)\hat g(k). \end{aligned} $$

(2.131)

assuming that the individual Fourier transforms, $\hat f(k)$ and $\hat g(k)$ do exist.

In the context of scale-invariance, it is useful to note that:

$$\displaystyle \begin{aligned} \mathrm{F}\left[\,f(ax) \right] = \frac{1}{|a|}\hat f\left(\frac{k}{a}\right). \end{aligned} $$

(2.132)

Similarly, the Fourier transform of a power law (see Table 2.7) is also very useful. In particular, it can be shown that, the Fourier transform of f(x) = |x|^−a, with 0 < a < 1, is given by another power-law.⁶¹ Namely,

$$\displaystyle \begin{aligned} \mathrm{F}( |x|{}^{-a}) = \varGamma(1-a) k^{a-1}, \end{aligned} $$

(2.133)

being Γ(x) Euler’s gamma function. Furthermore, it can be shown that if a function f(x) ∼|x|^−a with 0 < a < 1 for x →∞, then its Fourier transform $\hat f(k) \sim k^{a-1}$ for k → 0.

Appendix 2: Numerical Generation of Series with Prescribed Statistics

One of the most commonly used method generate time series with a prescribed set of statistics is the inverse function method. It is based on passing a time series u whose values follow a uniform distribution⁶² in [0, 1] through an invertible function. That is, we generate a collection of values using:

$$\displaystyle \begin{aligned} x = F^{-1}(u). \end{aligned} $$

(2.134)

The series of values for x is then distributed according to the pdf:

$$\displaystyle \begin{aligned} p(x) = \frac{dF}{dx}. \end{aligned} $$

(2.135)

To prove this statement, let’s consider the cumulative distribution function associated with p(x) (Eq. 2.4) at x = a, that verifies:

$$\displaystyle \begin{aligned} \mathrm{cdf} (a) = P\left (X \leq a \right) = P\left (F^{-1}(U) \leq a \right) = P\left ( U \leq F(a) \right) = F(a). \end{aligned} $$

(2.136)

The fourth step follows from the fact that the cumulative function being an increasing function of its argument. The last step is due to the fact that the cumulative distribution function of the uniform distribution is simply its argument! If F(x) = cdf(x), then Eq.2.135 follows after invoking Eq. 2.5.

Generating time series with prescribed statistics is thus reduced to knowing the inverse of the cumulative distribution function of the desired distribution.⁶³ For instance, in the case of the exponential pdf, $E_\tau (x) = \tau ^{-1}\exp (-x/\tau )$, the sought inverse is,

$$\displaystyle \begin{aligned} \mathrm{cdf} (x) = 1 - \exp({-x/\tau }) \Longrightarrow \mathrm{cdf}^{-1}(x) = - \frac{\log(1- u) }{\tau}. \end{aligned} $$

(2.137)

Thus, one can generate a series distributed with an exponential pdf of mean value τ by iterating:

$$\displaystyle \begin{aligned} e_\tau = \frac{ \log(1 - u)}{-\tau}, \end{aligned} $$

(2.138)

using uniformly-distributed u ∈ [0, 1].

Regretfully, many important distributions do not have analytical expressions for their cdfs. Much less of their inverse! This is the case, for instance, of the Gamma pdf (Eq. 2.70). Another important example is the Lévy distributions, that only have an analytic expression for their characteristic function (Eq. 2.32). Thus, the inverse function method is useless to generate this kind of data. Luckily, a formula exists, based on combining two series of random numbers, u distributed uniformly in [−π/2, π/2], and v distributed exponentially with unit mean, that is able to generate a series of l values distributed according to prescribed Lévy statistics. For α ≠ 1 the formula reads [24, 25]:

$$\displaystyle \begin{aligned} l_{[\alpha, \lambda, \mu, \sigma]} = \mu + \frac{\sigma \sin\left(\alpha u +K_{\alpha, \lambda}\right) }{ \left[\cos\left(K_{\alpha,\lambda}\right) \cos{}(u)\right]^{1/\alpha}} \left[\frac{ \cos\left((1 - \alpha)u - K_{\alpha, \lambda}\right) }{v}\right]^{1/\alpha - 1} \end{aligned} $$

(2.139)

with $K_{\alpha , \lambda } = \tan ^{-1} \left ( \lambda \tan \left (\pi \alpha /2\right )\right )$. For the case α = 1, one should use instead [25]:

$$\displaystyle \begin{aligned} l_{[1, \lambda, \mu, \sigma]} =\mu + \frac{2\sigma}{\pi}\left\lbrace \left( \frac{\pi}{2} + \lambda u \right)\tan u - \lambda \left[ \log\left(\frac{v\cos u}{\sigma\left(\frac{\pi}{2} + \lambda u\right)}\right) \right] \right\rbrace. \end{aligned} $$

(2.140)

These formulas reduce to much simpler expressions in some cases. For instance, in the Gaussian case (α = 2, λ = 0), Eq. 2.139 reduces to:

$$\displaystyle \begin{aligned} g_{[\mu, \sigma]} = \mu + 2 \sigma\sqrt{v} \sin{}( u) \end{aligned} $$

(2.141)

whilst for the Cauchy distribution (α = 1, λ = 0), Eq. 2.140 reduces to:

$$\displaystyle \begin{aligned} c_{[\mu, \sigma]} = \mu + \sigma\tan{}(u). \end{aligned} $$

(2.142)

Problems

2.1 Survival and Cumulative Distribution Functions

Calculate explicitly the cumulative distribution function and the survival function of the Cauchy, Laplace and exponential pdfs.

2.2 Characteristic Function

Calculate explicitly the characteristic function of the Gaussian, exponential and the uniform pdfs.

2.3 Cumulants and Moments

Derive the relations between cumulants and moments up to order n = 5.

2.4 Gaussian Distribution

⁶⁴Show that for a Gaussian pdf, $N_{[0,\sigma ^2]}(x)$, all cumulants are zero except for c ₂ = σ ². Also, show that only the even moments are non-zero and equal to m _n  = σ ⁿ(n − 1)!!.

2.5 Log-Normal Distribution

Calculate the first four cumulants of the log-normal distribution defined by Eq. 2.45. Show that they are given by: $c_1 = \exp (\mu +\sigma ^2/2)$; $c_2 = (\exp (\sigma ^2) -1 )\exp (2\mu +\sigma ^2)$; $S = (\exp (\sigma ^2) + 2)\sqrt {\exp (\sigma ^2) - 1}$ and $K = \exp (4\sigma ^2) +2\exp (3\sigma ^2) + 3\exp (2\sigma ^2) - 6$.

2.6 Poisson Distribution

Derive the limit of the binomial distribution B(k, n;p) for n →∞ and p → 0, keeping np = λ, and prove that it is the Poisson discrete distribution (Eq. 2.67).

2.7 Ergodicity

Consider the process defined by $y(t) = \cos (\omega t + \phi )$, where ω (the frequency) is fixed, but with a different values of ϕ (the phase) for each realization of the process. Under what conditions are the temporal and ensemble views of the statistics of y equivalent? Or, in other words, when does the process behave ergodically?

2.8 Numerical Estimation of pdfs

Write a code that, given a set of experimental data, $\left \{x_i, ~~i= 1, 2, \cdots N \right \}$, calculates their pdf using each of the methods described in this chapter: the CBS, CBC and the survival/cumulative distribution function methods.

2.9 Maximum Likelihood Estimators (I)

Calculate the maximum likelihood estimators of the Exponential pdf (Eq. 2.59) and the bounded power-law pdf (Eq. 2.111) using the methodology described in Sect. 2.5.

2.10 Maximum Likelihood Estimators (II)

Prove that, in the case of the Gumbel distribution (Eq. 2.51), the determination of the maximum likelihood estimators for μ and σ requires to solve numerically the pair of coupled nonlinear equations [26]:

$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat \mu &\displaystyle =&\displaystyle -\hat \sigma \log\left[\frac{1}{N}\sum_{i=1}^N\exp(-X_i/ \hat \sigma) \right] \end{array} \end{aligned} $$

(2.143)

$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat \sigma &\displaystyle =&\displaystyle \frac{1}{N}\sum_{i=1}^N X_i - \frac{ \displaystyle\sum_{i=1}^N X_i \exp(-X_i/\hat \sigma)} {\displaystyle \sum_{i=1}^N \exp (-X_i/ \hat \sigma)} \end{array} \end{aligned} $$

(2.144)

Write a numerical code that solves this equation using a nonlinear Newton method [18].

2.11 Synthetic Data Generation

Write a code that implements Eq. 2.139 to generate Lévy distributed synthetic data with arbitrary parameters. Use it to generate at least three different sets of data, and assess the goodness of Eq. 2.139 by computing their pdf.

2.12 Running Sandpile: Sandpile Statistics

Write a code that, given a flip time-series generated by the running sandpile code built in Problem 1.5, is capable of extracting the size and duration of the avalanches contained in it, as well as the waiting-times between them. Use the code proposed in Problem 2.8 to obtain the pdfs of each of these quantities.

2.13 Advanced Problem: Minimum Chi-Square Parameter Estimation

Write a code that, given $\left \{x_i, ~~i= 1, 2, \cdots N \right \}$, calculates the characteristic function of its pdf, ϕ ^data(k _j ), with k _j being the collocation points in Fourier space. Then, estimate the parameters of the Lévy pdf that best reproduces the data by minimizing the target function, $\chi _k^2(\alpha ,\lambda ,\mu ,\sigma )$, that quantifies the difference between the characteristic function of the data and that of the Lévy pdf (Eq. 2.32):

$$\displaystyle \begin{aligned} \chi^2_k(\alpha,\lambda,\mu,\sigma) = \sum_{i=1}^N \frac{ \left|\phi^{\mathrm{data}}(k_j) - \phi^L_{[\alpha,\lambda,\mu,\sigma]}(k_j) \right|{}^2}{\left| \phi^L_{[\alpha,\lambda,\mu,\sigma]}(k_j) \right|}, \end{aligned} $$

(2.145)

To minimize $\chi ^2_k$, the reader should look for any open-source subroutine that can do local optimization. We recommend using, for example, the Levenberg-Marquardt algorithm [18].