Statistics: descriptive measures and hypothesis testing

1. Population, sample and variables

This chapter is the one the deck describes, in its subtitle, as just what I'll need (not for statisticians!) — the statistics course embedded in an analytics course. It is built in two movements: descriptive statistics, quantities that summarise a data set, and inferential statistics, techniques that generalise from a sample to a population. The bridge between the two is the vocabulary that follows.

A population is the entire set of data, individuals, objects or results of interest. It is often too large to be fully analyzed, and it can be real or hypothetical — for example, the results of an experiment repeated countless times. A sample is a subset of the population. A sample can be random — each member is equally likely to be extracted — or a reasoned choice (non-probabilistic). Random selection seeks to ensure that the sample is representative of the population.

Variables are the quantities measured in a sample, and the first classification any analysis needs is their type:

Two further terms close the vocabulary. Descriptive statistics are quantities and techniques used to describe the characteristics of a data set, e.g. mean, standard deviation, box-plot. Parameters are quantities that describe the characteristics of a population; they are usually not known, and we want to make a statistical inference about them. Inferential statistics are techniques to analyze samples and generalize them to the population. And the gap between the two worlds has a name: the sample error is the difference between sample statistics and the values of the corresponding population parameters.

2. Frequency distributions, and the traffic counts

An (empirical) Frequency Distribution or Histogram for a continuous variable has a count of observations, grouped into predefined classes or groups. A Relative Frequency Distribution has the corresponding proportions of observations within the classes. A bar chart presents frequencies for a categorical variable — the point where the taxonomy of section 1 becomes practical.

The running example of the whole descriptive half of the deck is a traffic dataset: daily counts of vehicles on a motorway section, six months' worth of columns. The first column, ago1, holds 31 observations, one per day of August:

5012  4948  5077  4960  5010  4170  2916  4442  4301  4931  4533  4438  3707
2794  3148  4086  4355  4272  4389  3546  2827  4505  4598  4646  4406  4586
3801  2845  4919  4663  5167

import numpy as np
import pandas as pd
from scipy import stats            # to be used later
import matplotlib.pyplot as plt

df = pd.read_csv('traffico16.csv')   # dataframe (series)
npa = df['ago1'].to_numpy()          # numpy array
plt.hist(npa, bins=10, color='#00AA00', edgecolor='black')
plt.title(df.columns[0]); plt.xlabel('num'); plt.ylabel('days')
plt.show()

Grouped into ten classes of equal width, the counts give the relative frequency table. The two right-hand columns are exactly what the deck prints; the count column is reproduced here recomputed from the data, because the slide's printed version garbles the last three bins:

res = stats.relfreq(npa, numbins=10)
print(res[0])

3. Measures of central tendency: mean, median, mode

The central tendency measures indicate in which area of the range of permissible values the data are located. The common measures are three: the arithmetic mean, the median and the mode.

Let x₁, x₂, x₃, …, x_n be the measured values of a random variable X from a sample of cardinality n. The arithmetic mean is defined as

x̄ = (xᵢ + xᵣ + … + xₙ) / n          In Excel: MEDIA(dati)   In Python: np.mean(npa)

Example — some of the speeds measured on a motorway section are 151, 124, 132, 170, 146, 124, 113 km/h. The average is (151 + 124 + 132 + 170 + 146 + 124 + 113) / 7 = 137.14.

The median: n sample data are organized by increasing values, then the median is the middle value if n is odd, the average between the two middle values if n is even. The mode is the most frequently detected value.

median:  In Excel: MEDIANA(dati)   In Python: np.median(npa)
mode:    In Excel: MODA(dati)      In Python: stats.mode(npa)

Example, n odd — the speeds seen before, ordered, are 113, 124, 124, 132, 146, 151, 170. The median is the middle value: 132. Two travellers were driving at 124 km/h, so the mode is 124.

Example, n even — hotel quotes for a skiing holiday: 366, 327, 274, 292, 274, 230. Rearranged: 230, 274, 274, 292, 327, 366. The median is halfway between the two central values: (274 + 292) / 2 = 283. Two hotels asked the same amount: the mode is 274.

The deck then states the rule that motivates every robust statistic in the course: if the sample contains very high or very low values, the average tends to be distorted; the median is not affected by very large (or very small) values, so it is a better measure of centrality when the distribution is distorted. And the symmetry criterion closes the trio: if average = median = mode then the data are called symmetric.

4. Quartiles, percentiles and the IQR

Using the same principle of increasing ordering of data and their location, it is possible to define various quantiles — dividing into 4 intervals gives quartiles, dividing into 100 ranges gives percentiles. The 75th percentile, for example, is the value such that, in ascending ordering: 75% of the data has a lower value (stays on the left in the sort) and 25% has a higher value (stays on the right). Note: the median is the 2nd quartile and the 50th percentile.

The median divides a distribution into two halves. The first and third quartiles, Q1 and Q3, are defined symmetrically: 25% of the data is below Q1 (and 75% above), 25% is above Q3 (and 75% below). The inter-quartile range is the difference between the first and third quartile:

IQR = Q3 - Q1

Example — the ordered speeds 113, 124, 124, 132, 146, 151, 170 have Q1 = 124 (two values below it: 113, 124) and Q3 = 151 (two values above it: 151, 170). The Inter Quartile Range is 151 − 124 = 27.

In Excel: manual (difference between quartiles)
In Python: df['ago1'].quantile([0.25, 0.5, 0.75])   stats.iqr(npa)

The IQR matters beyond description: it is the scale on which the outlier fences of section 6 are built, and together with the median it survives skewness that would corrupt the mean and the standard deviation.

5. Measures of dispersion: range, variance, CV

Dispersion measures characterize how much the sample is scattered around the mean value, how variable the data are. The commonly used measures are: range; variance and standard deviation; coefficient of variation; inter-quartile range (seen before).

The Range of the sample is the difference between the largest and smallest value. It is easy to compute — e.g. for speeds min = 25, max = 203, so range = 178 km/h — and it is useful for defining scenarios, the best or the worst. But it is very sensitive to extreme values: one measurement error rewrites the whole measure.

In Excel: MAX(dati)-MIN(dati)   (RANGE in English means something else)
In Python: max(npa)-min(npa)     (range() is something else)

Variance, σ², is the arithmetic mean of the square of deviations from the mean:

σ² = ∑(xᵢ - x̄)² / n        (population mean known)

with a note that will matter for every exam-project comparison: the estimator is correct if the population average is known; otherwise (only the sample is known) it is better to divide by n − 1 and not by n.

In Excel: VAR.P(dati), VAR.C(dati)      In Python: np.var(npa)

The standard deviation σ is the square root of variance, and its advantage is that it has the same unit of measurement as the original variable x.

In Excel: DEV.ST.P(dati), DEV.ST.C(dati)    In Python: np.std(npa)

Example — the speed data worked by hand:

x̄ = 137.14, and dividing the sum of squared deviations by 7 gives σ = √(2304.86 / 7) = 18.14 km/h.

The coefficient of variation (CV), or relative standard deviation (RSD), is the standard deviation expressed as a percentage of the mean: CV = (σ / x̄) × 100%. It is not affected by multiplicative variations of scale, so it is useful when you want to compare distributions of variables measured on different scales. Example: CV = 100 × 19.6 / 137.1 = 14.3%.

In Python: stats.variation(npa)

The deck also shows the one-liner that produces the entire descriptive summary at once:

print(df['ago1'].describe())
#   count      31.000000
#   mean     4258.000000
#   std       720.147901
#   min      2794.000000
#   25%      3943.500000
#   50%      4438.000000
#   75%      4791.000000
#   max      5167.000000
#   Name: ago1, dtype: float64

Note that std here is the sample standard deviation (n − 1), the same value the widget of section 2 computes from the embedded data.

6. Box-plots, outliers and scatter-plots

A box-plot is a visual representation of a distribution based on: minimum, Q1, median, Q3, maximum — the five-number summary. It is useful for comparing large data sets side by side.

Example — speeds on rural roads: 62, 64, 68, 70, 70, 74, 74, 76, 76, 78, 78, 80. Here Q1 = (68 + 70) / 2 = 69, Q3 = (76 + 78) / 2 = 77, so IQR = 77 − 69 = 8. In Python: df.boxplot(column=['ago1', 'ago2']) or plt.boxplot(npa).

An outlier is an observation with very different value from those of the other data. It may be due to a measurement problem or may be indicative of a sub/population with abnormally high or low values. To represent them in a box-plot, redefine the lower and upper limits of the lines as:

Lower limit = Q1 - 1.5 × IQR
Upper limit = Q3 + 1.5 × IQR

Data may not reach these values; if there are data below the lower or above the upper limit, they are considered outliers. Outliers can be disturbing: they distort the average and increase variability. Elimination rules of thumb: in a "normal" sample a value should always be within 3 standard deviations of the mean, and often external values at 1.5–2 SD are discarded a priori.

Finally, the scatter-plot represents the relationship between two continuous variables; it is useful in the early stages of an investigation, to determine if there can be high correlation between the two, and it makes outliers evident. Python: plt.scatter(df['ago1'].sort_values(), df['set1'].sort_values()).

7. Distributions, PDF, CDF and the Gaussian

Measures of center and spread only tell a part of the story. We often need to look at the shape of the distribution by creating histograms that show the frequency of values. A sample of data will then be described by a distribution, i.e. by a function that describes the relationship between observations in the sample space; a distribution is a theoretical curve that models a histogram's shape, and we can use these curves to calculate estimated probabilities on which to base our conclusions. The function will fit the data with a modification of its parameters, such as mean and standard deviation.

Distributions are often described in terms of their density functions, which describe how the proportion of data or likelihood of the proportion of observations change over the range of the distribution. There are two kinds:

f(x) = (1 / (σ√(2π))) · exp( -x² / (2σ²) )     for mean 0; in general, (x - μ)² replaces x²

from scipy.stats import norm
dom = np.arange(-5, 5, 0.001)
mean = 0.0
plt.plot(dom, norm.pdf(dom, mean, 1),   label="std=1")
plt.plot(dom, norm.pdf(dom, mean, 0.5), label="std=0.5")
plt.plot(dom, norm.pdf(dom, mean, 2),   label="std=2")
plt.legend(); plt.show()
# the same calls with norm.cdf(dom, mean, std) plot the cumulative curves

When the generating distribution is unknown, the deck closes with two fitting tools that try every candidate and report the best by error: Fitter fits all or selected distributions and ranks them by SSE, and distfit returns the best distribution with its location/scale parameters and overlays the PDF on the data histogram.

from scipy import stats
from fitter import Fitter
data = stats.gamma.rvs(2, loc=1.5, scale=2, size=10000)   # empirical data (e.g. gamma)
f = Fitter(data)               # or specify: distributions=['gamma', 'normal', 'expon']
f.fit(); f.summary()           # top fits by SSE
print(f.fitted_param)          # dict of best params per distribution

# or
from distfit import distfit
data = np.random.normal(0, 2, 10000)
dfit = distfit()
results = dfit.fit_transform(data)
print(results)                 # best: e.g. 'norm' with loc/scale
dfit.plot()                    # PDF overlay

Which distribution family wins on your data decides which tests section 12 may use — parametric ones assume a shape, non-parametric ones do not.

8. Inference: sampling, the Gaussian and z-scores

The second deck opens with the problem statement of the entire part: how to collect only a limited number of data, a sample, and through their analysis reach general conclusions which can be extended to the whole population. To reach these conclusions, inference must be used: the ability to draw general conclusions (about the population or universe) using only a limited number of variable data (the sample).

Then comes the result that underwrites everything that follows. If a sampling operation is repeated 20 times from the same population, each time with a different random sample, 20 different averages will be obtained. Key result: all these sample averages tend to assume a normal distribution, even if the population of origin is not normally distributed. The random sampling process itself is a phenomenon that is normally distributed.

A normal distribution in a variable X with mean μ and variance σ² is a statistical distribution with a known probability function, defined on the domain x ∈ (−∞, ∞) — the formula is the one in section 7, and its shape is the bell curve. The deck then fixes the empirical percentages that make the curve usable:

The z-score (standard score, normal score) is a way of transforming every single value of a normal distribution into its standardized equivalent, specifying how many standard deviations the value is far from the population average: z = (x − μ) / σ. And two properties of probability on the curve are worth stating because the whole of hypothesis testing leans on them: the area under the curve represents the set of all possible cases, i.e. the total probability; probabilities are never referred to a point, but to an interval, and represent the ratio of all cases within that range to the total number of cases.

9. Confidence intervals

Confidence intervals (IC) provide a range of values within which we believe, with a certain level of confidence, that the true value falls — the unknown population parameter, not the sample statistic. For population averages the deck gives the two workhorse intervals:

95%   x̄ ± 1.96 · σ/√n
99%   x̄ ± 2.58 · σ/√n

10. Hypothesis testing: H₀ and the p-value

We have experimental data and we want to describe them in aggregate with descriptive statistics and distributions. It is assumed that the sample was extracted from a given population; we wonder how plausible this hypothesis is. The deck's examples are the questions that will be asked in the exam project: "I have a machine that launches a needsAttention alarm every two days, do I have to worry?" and "I ran three different learning algorithms on my data, are the results equivalent or is one better than the others? Are they doing better than what I am doing by hand?"

Statistical hypothesis testing is a logical-mathematical process that leads to the conclusion that the hypothesis of randomness cannot be rejected or that it can be rejected, through the calculation of probabilities of committing an error with these statements. The hypothesis that the result obtained with the experimental data is due only to chance is called the null hypothesis and is denoted by H₀; it states that the differences between two or more groups are essentially attributable to chance.

The question a test answers, in the deck's precise phrasing: assuming that the differences between groups of empirical observations are due to exclusively random factors, what is the probability that among all the possible alternatives the situation described arises from the data collected (or an even more extreme one)? The decision rule follows: if this probability is relatively high, conventionally equal to or greater than 5%, the differences are attributed to purely random factors (acceptance of the null hypothesis); if it is low, say lower than 5%, it is accepted as likely that the differences are due to non-random factors (non-acceptance of the null hypothesis).

Example — flip a coin 10 times consecutively. The deck tabulates the chances of each number of heads, and marks the verdicts:

The p column is the probability of exactly that outcome under H₀ (a fair coin): C(10, k)/2¹⁰, in percent. Eight heads in ten — 4.39% — lands below the 5% line and the coin becomes suspicious; seven heads — 11.72% — does not.

The p-value is then defined: p = probability of observing a value more extreme than that considered, if the null hypothesis is true. The lower the value of p, the less plausible H₀ becomes as an explanation of the data. The deck illustrates it with a histogram of AGE values framed by two green lines: the probability of having a value outside the green line range if the null hypothesis is true is < 5% — results outside the green lines have p < 0.05, results inside have p > 0.05.

The interval view and the significance view are the same coin: if the value of the null hypothesis lies within the 95% interval, the null hypothesis is accepted (p > 0.05); if it lies outside the 95% interval, the null hypothesis is rejected (p < 0.05).

11. Degrees of freedom, tails, and Type I and II errors

Two parameters characterise a test before any data is seen. Degrees of freedom: the number of points, elements or other units in the input data which are free to vary — often equal to the number of observations minus 1. Tails: tests on one tail are used for already oriented hypotheses; two-tailed testing is used in all other cases. The two-sided convention is the default, and it is why section 9's critical values come in symmetric pairs.

Then the two errors that every test can make, and their notation:

The complementary situations are named too: not rejecting a true H₀ has probability (1 − α), and rejecting a false H₀ has probability (1 − β) — the power of the test. The deck's illustration is medical: a temperature threshold to decide healthy versus sick. Type 1 error is declaring a healthy person sick (a false alarm, in alarm-management terms); Type 2 error is declaring a sick person healthy (a missed alarm). The two distributions overlap, so no threshold can avoid both.

Then a subtlety that explains why the whole course reports p-values and never β: if the null hypothesis is false then some other hypothesis, H₁, must be true. If we cannot specify this alternative hypothesis, it is not possible to determine the probability of committing a type II error. Often no single alternative can be identified, so only the level of significance α is considered, without fixing β: it is considered more appropriate to guard against the most serious error, type I. The assessment is motivated by the asymmetry of the two conclusions: the rejection of H₀ implies that the alternative hypothesis is true, while its acceptance implies that "there are not enough elements to reject it" — a test never proves H₀; it fails to disprove it.

12. Parametric vs non-parametric, and the test flowchart

The main distinction between tests is parametric versus non-parametric. Parametric tests are based on assumptions on the distribution of population parameters; a normal (Gaussian) distribution is usually assumed. They are more powerful — more likely to detect a statistically significant effect if one exists — but they can be misleading if the underlying assumptions are not satisfied. Non-parametric tests make no assumptions about the distribution of the population (they are also called distribution-free tests); they are usually based on the ranks of observations, i.e. their order number rather than the observations themselves. They have less power and are less flexible than parametric tests.

Non-parametric tests are justified in three situations the deck enumerates: 1) the variables have obvious deviations from the normal (strongly asymmetric or with more than one peak); 2) the sample is too small to understand if there is a normal distribution of data; 3) the observations are represented by ordinal rankings (e.g. severity of a disease from 1 to 4).

The tests themselves fall into two families: 1) relationships between variables — correlation, regression, chi square; 2) differences between variables — tests on difference in effectiveness of different solution approaches, e.g. t-test, Analysis of Variance (ANOVA), Wilcoxon. The second family is the one the exam project lives in, and its decision about paired or unpaired observations is decisive: if data was collected from the same subjects for each model, tests are "paired" — in ML, the data for all models are the same. Data from two independent groups is "unpaired", as in a control group versus a treatment group.

The deck's flowchart resolves which test to use. Compressed to a table:

13. Normality tests

Because parametric tests assume the Gaussian, the deck devotes a full sequence to checking it. Normality tests check whether a sample of data is distributed normally or not, and there are two ways to do it: graphical methods (simple, rude, unreliable) and analytical methods (more reliable tests, such as the Kolmogorov-Smirnov, Anderson-Darling or Shapiro-Wilk tests).

Method 1, graphic — distribution symmetry analysis (requires a sufficiently large sample, at least 50 points): 1. calculate mean (A), median (M), range (R) and standard deviation (σ) of the sample; 2. A and M must be close together, within 1% of R (a normal distribution is symmetrical, A = M); 3. apply the 68-95-99.7 rule: 68% of the data within σ of the mean, 95% within 2σ, 99.7% within 3σ; 4. if both checks pass, the sample distribution may be normal.

Method 2, graph — histogram comparison: compare the data histogram with the normal curve with the same mean and standard deviation. It must have a single maximum, be symmetrical, with degrading tails. (The deck walks through doing this in Excel via the "Data analysis" add-in; in Python:)

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(20)
data = 20 * np.random.randn(100) + 100
nbins = 10
plt.hist(data, nbins - 1)
plt.show()      # compare with a normal of the same mean and std

Method 3 — the Q-Q plot. Quantile-quantile plot: compare the sample data with those that would be expected from the distribution of interest (here the Gaussian). Each sample datum is coupled with a datum of the reference distribution (cumulative); the points are plotted in a scatter plot with the idealized value on the x-axis and the sample data on the y-axis. A perfect match for the distribution shows as a line of points at an angle of 45 degrees from the lower left to the upper right; point deviations from the line show a deviation from the expected distribution.

The theory behind the coupling: normally distributed data has the same area between two successive data — for 7 points, the area under the curve between two neighbouring points is 1/7 of the total. The deck's worked Excel example uses the seven values −4, −3, 0.8, 1.8, 3.9, 6.2, 6.5 (n = 7, mean = 1.74, standard deviation = 3.85): a "cum z-score" column places the points at equispaced cumulative probabilities centred on 0.5 and separated by 1/n, i.e. 0.07, 0.21, 0.36, 0.50, 0.64, 0.79, 0.93, and a "cum norm" column converts those probabilities back into data units with INV.NORM.N(prob, sample mean, sample dev.st.):

Plotting DATA against cum norm, the points should approximate the 45-degree line. In Python the whole procedure is one call:

from statsmodels.graphics.gofplots import qqplot
data = np.array([-4, -3, 0.8, 1.8, 3.9, 6.2, 6.5])
qqplot(data, line='q')
plt.show()

Method 4 — statistical tests. Many tests exist, testing different hypotheses: Shapiro-Wilk (among the most powerful), Anderson-Darling, Chi-Square, Kolmogorov-Smirnov (a general test of distribution goodness of fit), D'Agostino K-squared (checks normality on the basis of skewness and kurtosis — skewness a measure of symmetry, kurtosis a measure of how much the tails are above or below normal ones), and more. The deck runs two of them on the same generated sample:

import numpy as np
from scipy.stats import shapiro
np.random.seed(20)
data = 20 * np.random.randn(100) + 100
alpha = 0.05
stat, p = shapiro(data)
print('Shapiro=%.3f, p=%.3f, alpha=%.3f' % (stat, p, alpha))
if p > alpha:
    print('Gaussian sample (no rejection of H0)')
else:
    print('NON-Gaussian sample (reject H0)')
#   Shapiro=0.983, p=0.241, alpha=0.050  → no rejection of H0

from scipy.stats import kstest
stat, p = kstest(data, 'norm')
print('Kolmogorov=%.3f, p=%.3f, alpha=%.3f' % (stat, p, alpha))
if p > alpha:
    print('Gaussian sample (no rejection of H0)')
else:
    print('NON-Gaussian sample (reject H0)')
#   Kolmogorov=1.000, p=0.000, alpha=0.050  → reject H0

14. Student's t, Diebold–Mariano and critical differences

Student's t distribution is a family of distributions varying with k = degrees of freedom. It is the distribution that replaces the standard normal when the population standard deviation is unknown and must be estimated from the sample — for a sample of size n the test statistic follows a t with n − 1 degrees of freedom, and the t curve is slightly fatter than the Gaussian (larger critical values) exactly because the estimation adds uncertainty. As n grows, t approaches the normal; section 9's z = 1.96 is the large-sample limit of the t critical value.

The last topic is the one the deck positions as the operational payoff of the whole chapter: Diebold-Mariano (DM, 1995). The setting: several plausible methods to forecast a time series; differences in forecast accuracy between the methods. Is this outcome due to pure chance, is the difference statistically significant or not? If there are just two plausible models there are model-free tests (they do not consider the models that generated the forecasts) that compare the forecasts, for example Diebold-Mariano, but more advanced tests exist. They all use a null hypothesis H₀ stating that the two forecasts have the same accuracy.

Skipping the math, the decision rule is a direct application of everything above: the test statistic DM is normally distributed; H₀ is rejected if the DM statistic falls outside the range [−z_α/2, z_α/2], i.e. if |DM| > z_α/2, where z_α/2 is the z-value from the standard normal distribution corresponding to half of the desired α level. Example with α = 0.05: since DM is two-tailed, 0.05 is split into 0.025 for the upper tail and 0.025 for the lower tail; the z-value for 0.025 is −1.96 and for 1 − 0.025 = 0.975 it is 1.96. H₀ is rejected if the DM statistic falls outside [−1.96, 1.96].

The deck's example compares an MLP (200 epochs) with an LSTM (100 epochs) on the Box-Jenkins dataset, forecasting the last 12 periods with a one-period lookahead per record:

from dm_test import dm_test      # https://github.com/johntwk/Diebold-Mariano-Test

MLP    = [417.662, 381.815, 423.952, 469.055, 471.502, 547.532, 639.609,
          596.040, 477.471, 445.010, 368.674, 414.231]
LSTM   = [416.999, 390.999, 418.999, 461.000, 471.999, 535.000, 622.000,
          606.000, 507.999, 461.000, 390.000, 431.999]
actual = [417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432]

rt = dm_test(actual, MLP, LSTM, h=1, crit="MSE")
print(rt)
#   dm_return(DM=2.933300449662608, p_value=0.013609843048263027)

The verdict is read on two scales, and both say the same thing: DM = 2.933 > 1.96, so the statistic falls outside [−1.96, 1.96] and H₀ is rejected; equivalently p = 0.0136 < 0.05, the probability that the accuracy difference is due to pure chance. The MLP and the LSTM do not have the same accuracy on this series — the difference is statistically significant, not a sampling artefact. The deck notes the error criterion (crit="MSE") and the horizon (h=1, one-period lookahead) as the two choices that define what "accuracy" means in the comparison.

Finally, when more than two models are in play, pairwise tests are gathered into a single picture: critical difference diagrams include in one diagram the results of several pairwise comparisons (via scikit-posthocs or the source library). The deck's diagram, based on Wilcoxon-Holm signed rank tests, shows the ranks in the ordering of 5 classes: classes 3 and 5 performed best, with no significant difference between them, while classes 1, 2 and 4 perform worse and are equivalent among them. The horizontal segments connecting groups that are not significantly different are the diagram's whole message: significance is a graph, not a ranking.

Check your understanding