Machine Learning & Signals Learning

5 Descriptive Statistics

5.1 Sample and Population Variance

Preliminaries

We assume a uni-variate random experiment that is described by a real-valued random variable \(X\) with

\begin{align*} \E [X] &=\mu \\ \Var [X] &= \sigma ^2. \end{align*}

Let

\begin{equation} \bx =\{x_1,\dots ,x_n\} \end{equation}

be \(n\) observations of \(X\), where \(n\) may be fixed in advance or chosen arbitrarily.

5.1.1 Definition

Sample mean Given observations \(x_1,\dots ,x_n\), the sample mean \(\bar {\bx }\) is given by

\begin{equation} \bar \bx = \frac 1n\sum _{i=1}^n x_i \end{equation}

Sample variance Sample variance is given by

\begin{align} s_{\text {unbiased}}^2 =\frac 1{n-1}\sum _{i=1}^n\bigl (x_i-\bar {x}\bigr )^2 \end{align} This formulas is unbiased.

Population variance The more intuitive (biased) population variance formula is

\begin{equation} s_{\text {biased}}^2 =\frac 1n\sum _{i=1}^n\bigl (x_i-\bar {x}\bigr )^2\\ \end{equation}

Standard Deviation The sample standard deviation (std) is simply the square-root of the (unbiased or biased) sample variance, \(s\).

5.1.2 Repeated Experiments Characterization

The repeated experiments approach is based on

• • \(k\) experiments

  • – In each experiment \(j=1,\ldots ,k\), a sample \(\bx _j\) of \(n\) observations is drawn from the known random variable \(X\).

Mean

The average of repeated experiments values is itself approximately \(\mu \), so on average the sample mean

\begin{equation} \bar {\mu }_k = \frac {1}{k}\sum _{j=1}^{k}\bar {\bx }_j \end{equation}

recovers the true mean

\begin{equation} \lim _{k\rightarrow \infty }\bar {\mu }_k \rightarrow \mu \end{equation}

This principle is illustrated in Fig. 5.1.

Fluctuations of \(\bar {\mu }_k\):

• • As you gather more observations (higher \(n\)), \(\bar \bx \) tends to stabilize: it fluctuates less, and its empirical distribution concentrates around the true center \(\mu \).

• • Moreover, for large \(n\), the variability of \(\bar \bx \) scales like \(\sigma /\sqrt {n}\) (or \(\Var [\bar x]\approx \dfrac {\sigma ^2}{n}\)), meaning the more data we collect, the tighter our estimate of the process’s center becomes.

Variance

In the case of biased variance, the mean of the \(s_1^2,\ldots ,s_k^2\) values

\begin{equation} \bar {s}_k^2 = \frac {1}{k}\sum _{j=1}^{k}s_j^2 \end{equation}

systematically underestimates the true population variance \(\sigma ^2\),

\begin{equation} \lim _{k\rightarrow \infty }\bar {s}_{k,\text {biased}}^2 \rightarrow \dfrac {n-1}{n}\sigma ^2 \end{equation}

This inherent difference is termed bias. However, for the unbiased variance,

\begin{equation} \lim _{k\rightarrow \infty }\bar {s}_{k,\text {unbiased}}^2 \rightarrow \sigma ^2 \end{equation}

This principle is illustrated in Fig. 5.2.

MSE

The reason that biased estimator is still useful, is the when the empirical mean-squared error (MSE) is of interest. MSE is defined by

\begin{align} \mathrm {\overline {MSE}}_{\text {unbiased}} &= \frac {1}{k}\sum _{j=1}^k \left (s_{\text {unbiased},j}^2 - \sigma ^2\right )^2 \rightarrow \frac {2}{n-1}\,\sigma ^4\\ \mathrm {\overline {MSE}}_{\text {biased}} &= \frac {1}{k}\sum _{j=1}^k \left (s_{\text {biased},j}^2 - \sigma ^2\right )^2 \rightarrow \frac {2n-1}{n^2}\,\sigma ^4 \end{align} and

\begin{equation} \mathrm {\overline {MSE}}_{\text {unbiased}} \geq \mathrm {\overline {MSE}}_{\text {biased}} \end{equation}

5.1.3 Bias

The bias of a sample statistic is the systematic difference between its long-run average over repeated experiments and the corresponding true value of the random variable. A statistic is called unbiased when this difference is zero, and biased otherwise.

For example, for the sample variance,

\begin{equation} \mathrm {bias}(s^2) = \lim _{k\to \infty }\bar {s}_k^2 - \sigma ^2. \end{equation}

Therefore, \(s_{\text {unbiased}}^2\) is unbiased while \(s_{\text {biased}}^2\) has bias \(-\sigma ^2/n\).

The application of biased and unbiased estimators:

• • Biased: in ML tasks (e.g., loss function), optimal (maximum-likelihood estimation) for Gaussian distribution

• • Unbiased: statistics.

Code implementation default varies:

• • Python uses biased expression, e.g. numpy.std.

• • Matlab uses unbiased expression, e.g. std.

5.2 Standard Error and Confidence Bounds

5.2.1 Standard Error

The fluctuation of \(\bar x\) noted above has a name and a closed form. The standard error is the standard deviation of an estimate, as opposed to the standard deviation of the data the estimate was computed from.

For the sample mean of \(n\) uncorrelated samples of variance \(\sigma ^2\), the variances of the \(n\) terms add, while the \(1/n\) in front of the sum is squared,

\begin{equation} \label {eq-se-variance} \Var [\bar x] = \frac {1}{n^2}\sum _{i=1}^{n}\Var [x_i] = \frac {n\sigma ^2}{n^2} = \frac {\sigma ^2}{n} , \end{equation}

so that

\begin{equation} \label {eq-standard-error} \mathrm {SE}(\bar x) = \sqrt {\Var [\bar x]} = \frac {\sigma }{\sqrt {n}} . \end{equation}

The square root appears because it is the variances that add, not the amplitudes.

Eq. (5.15) is read as a budget: the precision of the mean improves only with the square root of the sample size, so halving the standard error costs four times the data. The population \(\sigma \) is rarely known, and in practice it is replaced by the sample standard deviation \(s_x\),

\begin{equation} \label {eq-standard-error-estimated} \widehat {\mathrm {SE}}(\bar x) = \frac {s_x}{\sqrt {n}} . \end{equation}

Table 5.2 shows the rule on simulated data. Read as standard errors, that is square-rooted, the three empirical variances of \(\bar x\) are \(0.445\), \(0.221\) and \(0.100\), against \(\sigma /\sqrt {n} = 0.447\), \(0.224\) and \(0.100\). Quadrupling \(n\) from \(5\) to \(20\) halves the standard error, exactly as (5.15) requires.

5.2.2 Confidence Bounds

The standard error is a scale, not yet a statement about where the true value lies. A confidence bound turns it into an interval: a range computed from the data that covers the true value in a fraction \(1-\alpha \) of repeated experiments, where \(\alpha \) is the accepted risk of missing it. The usual choice is \(\alpha =0.05\), a \(95\%\) interval.

Mean Let the \(n\) observations be i.i.d. \(\mathcal {N}(\mu ,\sigma ^2)\), with both parameters unknown. That is the difficulty: a statement about \(\mu \) appears to need \(\sigma \), which is no more available than \(\mu \) is. The standardized sample mean (2.43) is the way around it. In \(\left (\bar x-\mu \right )/\left (s_x/\sqrt {n}\right )\sim t_{n-1}\), every quantity on the left except \(\mu \) itself is computed from the sample, while the distribution on the right involves neither \(\mu \) nor \(\sigma \). That distribution is therefore known before a single observation is collected, which is exactly what the next step needs.

Since \(t_{n-1}\) is symmetric about zero, leaving \(\alpha /2\) of the probability in each tail means staying between \(\pm t_{n-1,\,1-\alpha /2}\),

\begin{equation*} \Pr \left (-t_{n-1,\,1-\alpha /2} \le \frac {\bar x - \mu }{s_x/\sqrt {n}} \le t_{n-1,\,1-\alpha /2}\right ) = 1-\alpha . \end{equation*}

This is still a statement about the statistic, not about \(\mu \), but the inequality can be solved for \(\mu \). Multiply all three parts by \(s_x/\sqrt {n}\), subtract \(\bar x\), then reverse the signs, which reverses the inequalities and leaves \(\mu \) alone in the middle,

\begin{equation} \label {eq-ci-mean-prob} \Pr \left (\bar x - t_{n-1,\,1-\alpha /2}\frac {s_x}{\sqrt {n}} \le \mu \le \bar x + t_{n-1,\,1-\alpha /2}\frac {s_x}{\sqrt {n}}\right ) = 1-\alpha . \end{equation}

Both endpoints are computable from the sample alone, so (5.17) is the confidence bound written out. Compactly,

\begin{equation} \label {eq-ci-mean} \bar x \pm t_{n-1,\,1-\alpha /2}\,\frac {s_x}{\sqrt {n}} , \end{equation}

with \(t_{n-1,\,1-\alpha /2}\) the upper-tail quantile of Sec. 2.2.5.

The interval is symmetric about \(\bar x\), and the factor it is built from, \(s_x/\sqrt {n}\), is precisely the estimated standard error (5.16). It therefore narrows as \(1/\sqrt {n}\), at the same rate and for the same reason as the standard error itself. Its width is set by the quantile, and the quantile is drawn from \(t_{n-1}\) rather than from the Gaussian because \(\sigma \) was estimated instead of known. Sec. 2.2.5 shows how that penalty shrinks as \(n\) grows.

Variance The variance follows the same argument, with (2.39) in place of (2.43). That result reads \((n-1)s_x^2/\sigma ^2\sim \chi ^2_{n-1}\), and it qualifies for the same reason: the unknown \(\sigma ^2\) sits on the left, while the distribution on the right involves no unknown at all. The chi-square density is not symmetric, so here the two tails need quantiles of their own,

\begin{equation*} \Pr \left (\chi ^2_{n-1,\,\alpha /2} \le \frac {(n-1)s_x^2}{\sigma ^2} \le \chi ^2_{n-1,\,1-\alpha /2}\right ) = 1-\alpha . \end{equation*}

Solving for \(\sigma ^2\) means taking reciprocals, which again reverses the inequalities and so exchanges the two endpoints: the larger quantile ends up in the lower bound,

\begin{equation} \label {eq-ci-variance} \left [\frac {(n-1)s_x^2}{\chi ^2_{n-1,\,1-\alpha /2}},\ \frac {(n-1)s_x^2}{\chi ^2_{n-1,\,\alpha /2}}\right ] . \end{equation}

The result is not symmetric about \(s_x^2\). Right skew places the two quantiles at unequal distances from the center, so the upper endpoint is the further one. Taking square roots of both endpoints gives the interval for \(\sigma \).

The two results are worth reading side by side. Ten observations locate the mean to within \(\pm 0.93\), while the variance is pinned only to a factor of seven between the endpoints, \(5.63/0.80\approx 7\). This is the concrete form of what the MSE comparison of Sec. 5.1 says in the abstract: the same data that places the center reasonably well barely constrains the spread, and a variance quoted from a small sample carries far less information than its digits suggest.

5.3 Mode, Median, Skewness and Kurtosis

5.3.1 Median

Median is a value or quantity lying at the midpoint of a frequency distribution of observed values or quantities, such that there is an equal probability of falling above or below it.

Equivalently, for a sample \(x_1,\dots ,x_n\) with order statistics \(x_{(1)}\le \cdots \le x_{(n)}\), the sample median is

\[ \mathrm {median}( x) = \begin {cases} x_{(\frac {n+1}2)}, & n\text { odd},\\[6pt] \dfrac {x_{(\frac n2)} + x_{(\frac n2 + 1)}}2, & n\text { even}. \end {cases} \]

5.3.2 Mode

For a finite sample \(x_1,\dots ,x_n\), the sample mode is any value(s) that occur(s) most often among the observations.

Multimodality: A distribution is called multimodal if it has more than one local maximum (mode) in its probability density function or histogram. Common cases include:

• • Unimodal: a single peak (e.g., normal distribution).

• • Bimodal: two distinct peaks, often indicating that the data is a mixture of two subpopulations.

• • Multimodal: three or more peaks.

Multimodality suggests that the data may come from a mixture of different underlying processes or groups.

5.3.3 Skewness

Skewness: Skewness measures the asymmetry of a distribution about its mean. The sample skewness is defined as

\[ g_1 = \frac {1}{n}\sum _{i=1}^{n}\!\left (\frac {x_i - \bar {x}}{s_n}\right )^{\!3}. \]

• • \(g_1 = 0\): the distribution is symmetric (e.g., normal distribution).

• • \(g_1 > 0\): the distribution is right-skewed (positively skewed) — the right tail is longer, and the mean is typically greater than the median.

• • \(g_1 < 0\): the distribution is left-skewed (negatively skewed) — the left tail is longer, and the mean is typically less than the median.

An illustration of the relationship between mean, median, and mode under different skewness is shown in Fig. 5.3.

5.3.4 Kurtosis

Kurtosis: Kurtosis measures the weight of the tails of a distribution, relative to a normal distribution of the same spread. The sample kurtosis is defined as

\[ g_2 = \frac {1}{n}\sum _{i=1}^{n}\!\left (\frac {x_i - \bar {x}}{s_n}\right )^{\!4}, \]

one power above the skewness \(g_1\). A normal distribution gives \(g_2 = 3\), so it is common to subtract this reference and report the excess kurtosis \(g_2 - 3\):

• • \(g_2 - 3 = 0\): mesokurtic — the tails match those of a normal distribution.

• • \(g_2 - 3 > 0\): leptokurtic — heavier tails, so observations far from the mean are more frequent than a normal distribution of the same \(\sigma \) predicts.

• • \(g_2 - 3 < 0\): platykurtic — lighter tails, so the extremes are effectively bounded.

An illustration of the three cases, at equal mean and equal \(\sigma \), is shown in Fig. 5.4.

Each deviation is divided by \(s_n\) before it is raised to the fourth power, so \(g_2\) carries no information about how wide the data are; it reports only how much of the data sits far from the center. This is what makes it worth quoting next to the std rather than instead of it. Two data sets with the same \(\bar {x}\) and the same \(s\) can still produce extreme values at very different rates: for the two distributions of Fig. 5.4a and Fig. 5.4b, a value beyond \(\bar {x}\pm 3s\) occurs with probability \(0.014\) and \(0.0027\) respectively, a factor of five. A limit placed at \(3s\) on the assumption of a normal distribution is therefore exceeded five times as often as intended when the data are leptokurtic.

Code implementation default varies here as well:

• • Python returns the excess form, e.g. scipy.stats.kurtosis with its default fisher=True.

• • Matlab returns the raw \(g_2\), e.g. kurtosis.

5.4 Monte-Carlo Simulation

The repeated experiments of Sec. 5.1 are an instance of a general numerical method, known as Monte-Carlo simulation. Instead of deriving the behavior of a statistic analytically, the random experiment is reproduced \(k\) times on a computer, and the resulting \(k\) values are summarized by the very same descriptive tools of this chapter and the next.

Formally, let \(g\) be any statistic computed from a sample of \(n\) observations, e.g. \(\bar {\bx }\), \(s^2_{\text {biased}}\) or the MSE. Its expected value is approximated by the average over \(k\) independent repetitions,

\begin{equation} \label {eq-monte-carlo} \E [g]\approx \bar {g}_k = \frac 1k\sum _{j=1}^{k}g_j, \end{equation}

where the law of large numbers guarantees \(\bar {g}_k \rightarrow \E [g]\) as \(k\rightarrow \infty \). The spread of \(g\) is estimated from the very same \(k\) values,

\begin{equation} \label {eq-mc-var} \Var [g]\approx s_g^2 = \frac 1k\sum _{j=1}^{k}\left (g_j-\bar {g}_k\right )^2,\qquad s_g=\sqrt {s_g^2}, \end{equation}

and \(\bar {g}_k-\theta \) estimates the bias with respect to a known true value \(\theta \) (Sec. 5.1.3).

The value of \(s_g\) is set by the sample size \(n\) and by the distribution of \(X\): it reports how much the statistic varies from one experiment to the next, and no number of repetitions reduces it. Making the statistic itself less variable requires a larger \(n\). The number of repetitions \(k\) controls something different, namely how accurately \(\bar {g}_k\) and \(s_g\) estimate the quantities they stand for.

The strength of the method is that (5.20) requires no derivation: it holds for any statistic, any distribution of \(X\), and any sample size \(n\), including all the cases where an analytical expression is unavailable or valid only for large \(n\).

The whole recipe is summarized in Fig. 5.5: the model feeds \(k\) independent experiments of \(n\) observations each, every experiment is reduced to a single value \(g_j\), and only the collection of the \(k\) values is interpreted.

Table 5.2 is precisely the output of this recipe, applied to the sample mean, the two variance estimators and their MSEs, with \(n\in \{5,20,100\}\) and \(k=5000\). Each of its entries condenses \(k\) stored values into a single number; Fig. 5.6 shows the values themselves behind one entry.

5.4.1 Number of Experiments

The simulated result is itself a random quantity, so it carries an error of its own, and this error is controlled by \(k\) alone.

Simulated mean The std of a simulated mean is

\begin{equation} \label {eq-mc-error-mean} \sqrt {\Var [\bar {g}_k]} = \frac {\sqrt {\Var [g]}}{\sqrt {k}}, \end{equation}

where \(\Var [g]\) is the spread (5.21) of the individual \(g_j\): the accuracy of a simulated mean is the spread of the statistic divided by \(\sqrt {k}\).

Simulated std and variance Similarly, the relative error of a simulated std \(s_g\) is

\begin{equation} \label {eq-mc-error-std} \frac {\sqrt {\Var [s_g]}}{s_g}\approx \frac {1}{\sqrt {2k}}, \end{equation}

while the relative error of a simulated variance \(s_g^2\) is twice as large,

\begin{equation} \label {eq-mc-error-var} \frac {\sqrt {\Var [s_g^2]}}{s_g^2}\approx \sqrt {\frac {2}{k}} = \frac {2}{\sqrt {2k}}. \end{equation}

The factor of two follows from squaring: a variance is the square of a std, and squaring doubles the relative fluctuation. A simulated variance is therefore always twice as uncertain as the simulated std obtained from the very same \(k\) values, and the two should not be quoted with the same accuracy.

All three decay as \(1/\sqrt {k}\), so halving the simulation error costs four times more experiments. Representative values are listed in Table 5.3 and illustrated in Fig. 5.7 for the simulated mean and in Fig. 5.8 for the simulated std.

The consequence is that a simulated value is never exact, and the number of correct digits it carries is set by \(k\). In Table 5.2, for instance, \(k=5000\) yields a variance error of about \(2\%\) (5.24), which is why the \(n=100\) row reads \(\bar s^2_{\mathrm {unbiased}}=1.01\) rather than the exact \(1.00\).

Derivation (*) Where do the three formulas come from? Each one answers the same question: if the whole simulation is repeated with a new seed, how far does the reported number move?

Start with the mean. The stored values \(g_1,\ldots ,g_k\) are independent, so their variances add up, and dividing their sum by \(k\) divides the variance by \(k^2\),

\begin{equation*} \Var [\bar {g}_k] = \frac {1}{k^2}\sum _{j=1}^{k}\Var [g_j] = \frac {\Var [g]}{k}. \end{equation*}

The square root of this is (5.22). It is the familiar \(\sigma ^2/n\) rule of the sample mean (Sec. 5.2.1), applied to \(k\) experiments instead of \(n\) observations. Only independence was used, so this rate is valid for any distribution of \(g\). This is also the one rate given as an absolute value instead of a percentage: a percentage would divide by \(\E [g]\), which is zero whenever the simulation checks a bias against zero.

The variance needs one known result. For Gaussian values, the sum of the squared deviations follows a chi-square distribution (Sec. 2.2.4), whose variance is twice its number of degrees of freedom. With \(k\) values this gives

\begin{equation*} \Var [s_g^2] \approx \frac {2\sigma _g^4}{k}, \end{equation*}

and dividing by \(s_g^2\approx \sigma _g^2\) turns it into the percentage (5.24).

The std needs no new assumption, only a square root. A small relative change of the variance shrinks to half of itself after the square root is taken,

\begin{equation*} s_g = \sigma _g\sqrt {1+\varepsilon }\approx \sigma _g\left (1+\frac {\varepsilon }{2}\right ), \end{equation*}

so a variance known to \(\pm 10\%\) gives a std known to \(\pm 5\%\). Halving \(\sqrt {2/k}\) gives \(1/\sqrt {2k}\), which is (5.23).

Note that only the middle step assumed a Gaussian. When \(g\) has heavier tails, rare large values inflate the squared deviations, and both the variance and the std are then less accurate than these formulas promise. The mean is unaffected. Whether the variance is normalized by \(k\) or by \(k-1\) makes no visible difference at these accuracies.

5.4.2 Reproducibility

A computer draws pseudo-random numbers, produced by a deterministic algorithm from an initial value called the seed. Fixing the seed makes the whole simulation exactly repeatable, which is a requirement for a reported result, not an option.

• • Python: numpy.random.default_rng(seed).

• • Matlab: rng(seed).

5.5 Exercises

5.1 (descriptive statistics of a small sample) A sensor is read eight times, giving

\[ 2,\ 4,\ 4,\ 4,\ 5,\ 5,\ 7,\ 9 . \]

• (a) Compute the sample mean \(\bar {x}\) and the sum of squared deviations \(\sum _i(x_i-\bar {x})^2\).

• (b) Compute the biased variance and its std, and the unbiased variance and its std.

• (c) Find the median, the mode and the range.

• (d) By what factor do the two variances differ? For which \(n\) does the choice between them change the reported value by no more than \(1\%\)?

• (e) Estimate the standard error of the mean, and give the \(95\%\) confidence interval for \(\bar {x}\). Compare the width of that interval with the range found in (c). Use \(t_{7,\,0.975} = 2.365\).

5.2 (simulation: both variance estimators) Verify by simulation the bias \(-\sigma ^2/n\) of Sec. 5.1.3, and the price paid for removing it.

Draw \(k=10^4\) independent experiments, each consisting of \(n=10\) samples from \(\mathcal {N}(0,1)\) (randn). For every experiment compute both variance estimators, and then across the \(k\) experiments report:

• (a) the mean of each estimator, compared with \(\sigma ^2\) and with \(\frac {n-1}{n}\sigma ^2\);

• (b) the MSE of each estimator about the true \(\sigma ^2\), compared with

  \[ \overline {\mathrm {MSE}}_{\text {unbiased}} = \frac {2}{n-1}\sigma ^4, \qquad \overline {\mathrm {MSE}}_{\text {biased}} = \frac {2n-1}{n^2}\sigma ^4 ; \]

• (c) which of the two estimators you would rather use, and why the answer is not the unbiased one.

Fix the seed so the run is reproducible, and lay the results out as in Table 5.2.