Statistics: descriptive measures and hypothesis testing

1. Population, sample and variables

Descriptive statistics summarize observed data; inferential statistics use a sampling design or model to make statements beyond those observations. The distinction matters even before choosing a plot.

A population is the set of units or outcomes the analysis concerns; a sample is the observed subset. In simple random sampling without replacement, all subsets of the chosen size are equally likely. Probability sampling more generally can use unequal, known inclusion probabilities. Non-probability selection does not provide that design-based foundation. Random selection alone does not ensure that a realized sample is representative or remove coverage, nonresponse and measurement problems. See Statistics Canada’s sampling guide.

A variable records a characteristic of an observational unit. Its type concerns the meaning of the values, not simply how they are stored:

A statistic is calculated from a sample; a parameter describes the population or model of interest. For a particular sample, sampling error is the difference between an estimator and the corresponding population quantity. It is not a synonym for all measurement errors, selection bias or forecast errors.

2. Frequency distributions and the traffic counts

A frequency distribution records counts for values, categories or numeric intervals; relative frequencies divide by the declared total. A histogram groups quantitative observations, including discrete counts, into ordered numeric bins. A categorical bar chart compares categories. For a histogram, specify the bin edges and endpoint convention: the bin count alone is insufficient. For unequal widths, use count/width or relative frequency/width if bar area is to represent frequency or probability. Changing bins does not change the underlying observations.

The deck’s ago1 column contains 31 observations, one daily vehicle count for August. These are discrete quantitative data, even though a histogram groups them into intervals:

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

A reproducible histogram requires both a bin count and edges. NumPy’s histogram documentation uses the observed minimum and maximum by default. SciPy’s relfreq documentation instead expands each end by (max − min) / [2(k − 1)] for k bins. These defaults need not produce the same counts.

Use the same edges in Python

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

data = np.array([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])
k = 10
pad = (data.max() - data.min()) / (2 * (k - 1))
edges = np.linspace(data.min() - pad, data.max() + pad, k + 1)
counts, _ = np.histogram(data, bins=edges)
relative = counts / data.size
cumulative = np.cumsum(counts) / data.size
print(counts)  # [4 1 0 2 1 2 8 5 6 2]
print(counts.sum())  # 31
res = stats.relfreq(data, numbins=k,
                    defaultreallimits=(edges[0], edges[-1]))
assert np.allclose(relative, res.frequency)
plt.hist(data, bins=edges, color='#1546B8', edgecolor='#F3EFE3')
plt.xlabel('Vehicles per day'); plt.ylabel('Days')
plt.show()
# For the min–max convention, use edges = np.linspace(data.min(), data.max(), k + 1).
# Recompute counts, proportions and the plot together after changing the edges.

Relative frequency is count / n; frequency density additionally divides by bin width. With unequal bin widths, plotting raw counts as heights makes the bar areas misleading. See the NIST histogram guide for normalization choices. Neither a histogram nor its binning establishes independence or a distributional model.

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

Mean, median and mode summarize different aspects of location. For the numerical samples here, the arithmetic mean is x̄ = Σxᵢ/n; the median is the middle sorted value for odd n and the average of the two middle values for even n. Modes are values attaining the greatest observed frequency; there can be ties. A density’s peak and a histogram’s tallest bin are not the same object as a sample mode. See NIST’s location guide.

The seven speeds are 151, 124, 132, 170, 146, 124, 113 km/h: their sum is 960, so x̄ = 960/7 = 137.142857…; the sorted values are 113, 124, 124, 132, 146, 151, 170. The six hotel quotes are 366, 327, 274, 292, 274, 230: the two central sorted values are 274 and 292, giving median 283.

Python: np.mean(x), np.median(x), and stats.mode(x). SciPy mode returns only one value when modes tie; use counts to retain all modes, as in the reproduction below. Italian Excel names include MEDIA, MEDIANA, and MODA.SNGL/MODA.MULT for single/multiple modes (the slides use legacy MODA); see the Microsoft function reference.

Robustness is not a claim that one summary is always better

Replacing one value x by x + Δ changes the sample mean by exactly Δ/n. The median depends on ranks and is less sensitive to moving an already extreme observation farther out, but it can change when observations cross the middle or enough values are replaced. Section 6 demonstrates this explicitly. Choose a summary for the question: expected total workload and typical daily workload need not have the same answer. Extreme but valid observations are not automatically errors.

For finite data, symmetry about c means that reflection x → 2c−x preserves values and multiplicities. For a distribution, the reflected random variable must have the same distribution; its mean need not exist. Equality of mean, median and mode is neither a definition nor a sufficient test of symmetry. A symmetric distribution or sample may also have modes away from its center.

import numpy as np

def transformed(a):
    u = (np.arange(121) - 60) / 60
    return 100 + 60*u if a == 0 else 100 + 60*np.expm1(a*u)/a

x = transformed(1.5)  # deterministic grid, not a random sample
edges = np.linspace(x.min(), x.max(), 15)
counts, _ = np.histogram(x, bins=edges)
print(x.mean(), np.median(x))
print(counts)
print(np.flatnonzero(counts == counts.max()) + 1)  # ALL highest-count bins
assert len(np.unique(x)) == 121  # no unique sample mode
# The central grid point u=0 maps to x=100 for every a.

4. Quartiles, percentiles and the IQR

Quantiles describe positions in an ordered distribution; percentiles use p/100, and Q1, Q2, Q3 correspond to p = 0.25, 0.5, 0.75. With the usual numerical median convention, Q2 is the median. Ties and finite samples prevent a guarantee that exactly 25% or 75% of observations are strictly below a reported quartile. Sample quantiles also depend on the selected estimator.

One population convention is Q(p) = inf{x: F(x) ≥ p} for 0 < p < 1. Interpolated sample estimators need not equal that inverse applied to the empirical CDF. NumPy’s linear method uses h = (n−1)p on the zero-based sorted array and interpolates between floor(h) and ceil(h). The slides instead use medians of halves in the worked examples; for seven values the middle observation is excluded from both halves.

IQR = Q3 − Q1. The slide’s 151 − 124 = 27 is correct for its convention; the linear result is 148.5 − 124 = 24.5. Only the value 113 is strictly below Q1 = 124; the two observations equal to 124 are not “below” it. Likewise 151 is not strictly above Q3 = 151. This is why strict percentage claims are misleading with ties and seven observations.

np.quantile(x, [.25, .5, .75], method='linear')
pd.Series(x).quantile([.25, .5, .75], interpolation='linear')
stats.iqr(x, interpolation='linear')  # 24.5 for these seven speeds

The IQR is a spread measure based on the central ranks, not squared distance from the mean. It is less sensitive to moving tail observations farther out, but it is not invariant under arbitrary data changes. Use the same quartile convention when calculating the fences in section 6.

5. Measures of dispersion: range, variance, CV

Spread has several definitions. Range = max − min uses the extremes; IQR uses quartiles; variance averages squared deviations from a specified center. Their sensitivity and units differ. The slide’s separate range example has min 25 and max 203 km/h, hence range 178; the seven motorway speeds instead have range 170 − 113 = 57 km/h. Python np.ptp(x) or max(x)-min(x) calculates a range; range() does not.

Separate descriptive spread from a variance estimator

Observed mean: x̄ = Σxᵢ / n
SS = Σ(xᵢ − x̄)²
Descriptive variance: vₙ = SS / n
Sample variance:      s² = SS / (n−1), n > 1
If the true mean μ is known: Σ(xᵢ − μ)² / n

vₙ describes the observed values; it does not require knowing a population mean. For an i.i.d. sample with finite population variance, s² is an unbiased estimator of that variance. If the true μ is known, the final formula is unbiased under that model, but it centers on μ, not on x̄. Do not label SS/n around the sample mean as the known-mean formula. These statements do not automatically apply to dependent traffic days or supply a valid confidence interval.

Standard deviation is the square root of the chosen variance. It has the original measurement units; variance has squared units. Taking the square root of an unbiased variance estimator does not generally produce an unbiased standard-deviation estimator.

SS = 16134/7 = 2304.857143…; calculations use the unrounded mean 960/7. Cells and totals are rounded independently for display. In particular, √(SS/7) = 18.145669…, which rounds to 18.15 km/h at two decimals; the slide prints 18.14.

NumPy var and std default to ddof=0 (divisor n); use ddof=1 for n−1. pandas Series.std defaults to ddof=1. Defaults belong to individual functions, not to a library as a whole. Italian Excel uses VAR.P/DEV.ST.P for the n convention and VAR.C/DEV.ST.C for n−1.

CV requires a meaningful scale and a declared denominator

For positive ratio-scale quantities with a meaningful zero and positive mean, CV percent = 100 × SD/x̄. Its relative-spread interpretation needs that scale: multiplying all observations by the same positive constant leaves it unchanged, but adding an offset generally changes it. A Celsius/Fahrenheit conversion is not a pure rescaling. CV is undefined at zero mean and unstable near zero; it is not a universally meaningful comparison of arbitrary variables. See the NIST CV guide.

SciPy stats.variation returns SD/mean, not a percentage, and defaults to ddof=0. It does not take the absolute value of the mean, so it can return a negative number without establishing a meaningful relative-spread interpretation. For the speeds, multiply by 100 to obtain 13.23% with ddof=0 or 14.29% with ddof=1. The slide’s 14.3% agrees with the latter after rounding; its plain stats.variation(npa) call would not reproduce that convention.

import numpy as np
import pandas as pd
from scipy import stats

x = np.array([151, 124, 132, 170, 146, 124, 113], dtype=float)
mean = x.mean()
ss = np.sum((x - mean)**2)
print("Mean, median:", mean, np.median(x))
values, counts = np.unique(x, return_counts=True)
print("All sample modes:", values[counts == counts.max()])
print("Linear quartiles:", np.quantile(x, [.25, .5, .75], method="linear"))
s = np.sort(x)
print("Median-of-halves quartiles:", np.median(s[:3]), np.median(s[4:]))
print("Sum of squared deviations:", ss)
for ddof in [0, 1]:
    print("ddof, variance, SD, CV percent:", ddof,
          np.var(x, ddof=ddof), np.std(x, ddof=ddof),
          100 * stats.variation(x, ddof=ddof))
assert np.isclose(stats.variation(x), np.std(x, ddof=0)/mean)
assert np.isclose(pd.Series(x).std(), np.std(x, ddof=1))
assert np.isclose(stats.iqr(x), 24.5)
# Same numeric array, no missing values; defaults are function-specific.

For the 31 August counts, pd.Series(data).describe() reports mean 4258, sample SD 720.147901, min 2794, linear quartiles 3943.5 / 4438 / 4791, and max 5167. These summarize the recorded days; the output does not prove an i.i.d. sampling design. The frequency and box widgets use this same data source.

6. Box-plots, outliers and paired scatter-plots

A box-plot summarizes ordered quantitative data using a box from Q1 to Q3 and a median line. Specify its convention: a min–max plot extends its whiskers to the extremes; the 1.5·IQR variant used here ends them at the most extreme observations still inside the fences. Whiskers are observed values; fences are calculated thresholds. Neither is a confidence interval.

IQR = Q3 − Q1
Lower fence = Q1 − 1.5 × IQR
Upper fence = Q3 + 1.5 × IQR
Flag x only if x < lower fence or x > upper fence.

The slide 30 data are 62, 64, 68, 70, 70, 74, 74, 76, 76, 78, 78, 80 km/h. Taking the median of each six-value half gives Q1 = (68 + 70)/2 = 69 and Q3 = (76 + 78)/2 = 77. Do not assume every library reproduces those quartiles: NumPy’s linear method interpolates at index (n−1)p and gives 69.5 and 76.5 for the same sample. Both conventions below have median 74, whiskers 62 and 80, and no flagged values.

import numpy as np
import matplotlib.pyplot as plt

speeds = np.array([62, 64, 68, 70, 70, 74, 74, 76, 76, 78, 78, 80])
s = np.sort(speeds)
# Reproduce this slide's median-of-halves convention explicitly.
q1, med, q3 = np.median(s[:6]), np.median(s), np.median(s[6:])
iqr = q3 - q1
inside = s[(s >= q1 - 1.5 * iqr) & (s <= q3 + 1.5 * iqr)]
fliers = s[(s < q1 - 1.5 * iqr) | (s > q3 + 1.5 * iqr)]
fig, ax = plt.subplots()
ax.bxp([dict(q1=q1, med=med, q3=q3,
             whislo=inside.min(), whishi=inside.max(), fliers=fliers)],
       orientation="horizontal")
ax.set_xlabel("Speed (km/h)")
print(q1, med, q3, iqr)  # 69.0 74.0 77.0 8.0
print(np.quantile(s, [.25, .5, .75], method="linear"))
# [69.5 74.  76.5] — a different, explicitly chosen convention
plt.show()

A flag is a prompt to investigate, not to delete

An extreme observation can be a recording error, a legitimate rare event, or evidence that the model or grouping needs attention. Check its provenance and context; correct established errors and document justified exclusions. Otherwise retain the observation and consider robust summaries or sensitivity analyses. NIST distinguishes flagging from accommodation and formal identification.

The slides’ “always within 3 SD” wording is not a guarantee: a normal distribution has unbounded tails. For a single observation under known normal parameters, P(|X−μ| > 3σ) ≈ 0.0027, not zero. Estimated sample z-scores are a different calculation. Neither 1.5·IQR nor a 1.5–2 SD cutoff licenses automatic removal; trimming valid tails changes the data distribution and can bias the analysis.

Scatter-plots must preserve the observational pairs

A scatter-plot places each matched pair (xᵢ, yᵢ) at its two numeric coordinates; discrete counts are allowed too. It can reveal nonlinear structure, clusters and unusual observations, not just linear correlation. Association does not establish causation, and zero Pearson correlation does not establish independence.

# Use this only if each existing row represents a meaningful matched unit.
# Otherwise first join on a verified observational key, not row position.
paired = df[['ago1', 'set1']].dropna()  # Remove incomplete PAIRS together.
plt.scatter(paired['ago1'].to_numpy(), paired['set1'].to_numpy())
# No independent sort_values(); sorting whole rows is harmless.

For counts from different months, matching equal day numbers is itself a design choice, not simultaneous observation. The available slide extraction does not establish a verified observational key for the two series. The illustration below therefore uses explicitly synthetic pairs, not an invented August–September correlation.

import numpy as np
import matplotlib.pyplot as plt

# Synthetic observations A–E; not the course's traffic data.
xy = np.array([[1, 4], [2, 1], [3, 3], [4, 5], [5, 2]])
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
axes[0].scatter(xy[:, 0], xy[:, 1])
wrong = np.column_stack((np.sort(xy[:, 0]), np.sort(xy[:, 1])))
axes[1].scatter(wrong[:, 0], wrong[:, 1])  # Deliberately WRONG pairing
axes[0].set_title("Original pairs")
axes[1].set_title("Sorted separately: artificial association")
for ax in axes:
    ax.set(xlim=(.5, 5.5), ylim=(.5, 5.5), xlabel="x", ylabel="y")
    ax.set_aspect("equal")
print(np.corrcoef(xy.T)[0, 1])     # 0.0
print(np.corrcoef(wrong.T)[0, 1])  # 1.0 (up to floating-point rounding)
plt.show()

7. Distributions, PDF, PMF, CDF and the Gaussian

A probability distribution assigns probabilities to events. An empirical distribution describes observed frequencies; a theoretical model proposes probabilities for a random variable X. A histogram summarizes a sample using chosen bins, but a distribution is not necessarily a smooth curve fitted to a histogram. Distinguish counts, probability masses, densities and cumulative probabilities.

f(x) = exp(−(x−μ)² / (2σ²)) / (σ√(2π)),  −∞ < x < ∞

One interval, two equivalent probability calculations

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

x = np.linspace(-5, 5, 1001)  # display window, not the full support
fig, axes = plt.subplots(1, 2)
for sigma in [0.5, 1, 2]:
    axes[0].plot(x, norm.pdf(x, loc=0, scale=sigma), label=f"SD={sigma}")
    axes[1].plot(x, norm.cdf(x, loc=0, scale=sigma), label=f"SD={sigma}")
for ax, name in zip(axes, ["Density f(x)", "Cumulative probability F(x)"]):
    ax.set(xlabel="x", ylabel=name)
    ax.legend()
fig.tight_layout()
plt.show()

# The paired teaching figure uses a narrower normal: SD=0.25.
rv = norm(loc=0, scale=0.25)
print("Density at zero:", rv.pdf(0))  # 1.595769..., not a probability
print("P(-0.25 < X <= 0.25):", rv.cdf(0.25) - rv.cdf(-0.25))
# For this continuous model, P(X=0)=0.

Fitting compares candidates; it does not certify a model

The slides introduce Fitter and distfit. These examples retain their two separate synthetic samples: gamma with shape 2, location 1.5 and scale 2; then normal with mean 0 and SD 2. Both are seeded for reproduction, use 100 density-histogram bins and compare only norm, gamma and expon. The SciPy candidate name is norm, not normal. Parameters are fitted freely; the generating parameters are not supplied to the fitting tools.

For the settings below, SSE/RSS is Σⱼ(hⱼ − f̂(cⱼ))², comparing histogram density hⱼ with the fitted PDF at bin center cⱼ. Lower is better for that sample, binning and candidate set. This is not a p-value or an integrated probability; changing the bins can change the ranking. Candidate failure is checked explicitly. Neither sample represents the traffic observations elsewhere in the chapter.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from fitter import Fitter

rng = np.random.default_rng(20260908)
data = stats.gamma.rvs(2, loc=1.5, scale=2, size=10000, random_state=rng)
candidates = ["norm", "gamma", "expon"]  # SciPy uses norm, not normal
f = Fitter(data, distributions=candidates, bins=100, density=True, verbose=False)
f.fit(max_workers=1, progress=False)
scores = f.df_errors["sumsquare_error"].reindex(candidates)
assert np.isfinite(scores.to_numpy()).all(), "At least one fit failed"
print(scores.sort_values())
print(f.get_best(method="sumsquare_error"))
name = scores.idxmin()
params = f.fitted_param[name]  # shape parameters (if any), loc, scale
grid = np.linspace(data.min(), data.max(), 501)
fig, ax = plt.subplots()
ax.hist(data, bins=100, density=True, alpha=0.25)
ax.plot(grid, getattr(stats, name).pdf(grid, *params), label=name)
ax.set(xlabel="Synthetic value", ylabel="Density")
ax.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from distfit import distfit

rng = np.random.default_rng(20260908)
data = rng.normal(loc=0, scale=2, size=10000)  # a separate synthetic sample
candidates = ["norm", "gamma", "expon"]
dfit = distfit(method="parametric", distr=candidates, stats="RSS", bins=100,
              mhist="numpy", n_boots=None, n_jobs=1, verbose="warning")
dfit.fit_transform(data)
scores = dfit.summary.set_index("name")["score"].reindex(candidates).astype(float)
assert np.isfinite(scores.to_numpy()).all(), "At least one fit failed"
print(scores.sort_values())
print(dfit.model["name"], dfit.model["params"])
grid = np.linspace(data.min(), data.max(), 501)
fig, ax = plt.subplots()
ax.hist(data, bins=100, density=True, alpha=0.25)
ax.plot(grid, dfit.model["model"].pdf(grid), label=dfit.model["name"])
ax.set(xlabel="Synthetic value", ylabel="Density")
ax.legend()
plt.show()

These snippets were executed with Fitter 1.8.0, distfit 2.0.2, NumPy 2.3.3, SciPy 1.16.2 and Matplotlib 3.10.7. The overlays deliberately show only the density fit, not automatically generated probability thresholds labelled as confidence intervals. In distfit, n_boots=None requests no bootstrap validation.

A winning candidate is not proof of the true family. In-sample fitting and model selection must be accounted for when calibrating goodness-of-fit tests; a standard known-parameter KS p-value is not automatically valid after estimating those parameters from the same sample. Parametric does not mean “always normal”, and nonparametric does not mean “assumption-free”. Choose a test for the question and sampling design, including pairing and dependence: see section 12 and the normality caveats in section 13. A density-fitting leaderboard alone does not authorize a test.

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

Inference estimates population/model quantities and quantifies uncertainty under a stated sampling design or probability model. A statistic such as x̄ varies from sample to sample; the target μ stays fixed within that model. The realized error x̄−μ is not the distance between arrows in a diagram. Random variation, systematic bias and measurement error are different issues.

Sample size n is not the number of repetitions R

In the classical i.i.d. CLT, X₁,…,Xₙ are independent and identically distributed with finite mean μ and finite positive variance σ². As the number of observations within each sample grows, the standardized mean converges in distribution to N(0,1):

E[X̄ₙ] = μ;   Var(X̄ₙ) = σ²/n;   SD(X̄ₙ) = σ/√n
Zₙ = (X̄ₙ − μ)/(σ/√n)  →  N(0,1) as n → ∞

Repeating the sampling R times helps visualize the distribution at fixed n. Increasing R does not make that distribution more normal. Twenty repetitions are not a CLT threshold, and no universal n guarantees an accurate approximation for every distribution and tail probability. Independent normal observations give an exactly normal mean for every n. See the MIT CLT and standardization notes.

σ/√n is a standard error: the SD of the estimator across samples, not the spread σ of individual observations. With dependence, Var(X̄) = [Σ Var(Xᵢ) + 2Σᵢ<ⱼ Cov(Xᵢ,Xⱼ)]/n². Covariances cannot simply be discarded. Sampling without replacement from a finite population needs its design-specific correction. Shared train/test splits, cross-validation and consecutive traffic days do not become independent merely by calling them repetitions.

Finite variance matters: an average of independent standard Cauchy observations is again standard Cauchy, not a concentrating normal mean; its population mean and variance do not exist. Other limit theorems have different conditions. This CLT does not license every test on arbitrarily dependent or heavy-tailed losses.

import numpy as np
from scipy.stats import gamma, norm

z = np.linspace(-4, 4, 801)
for n in [4, 16, 64]:
    xbar = 1 + z/np.sqrt(n)
    exact = gamma.pdf(xbar, a=n, scale=1/n)/np.sqrt(n)
    reference = norm.pdf(z)
    print(n, "SD of mean:", 1/np.sqrt(n), "skewness:", 2/np.sqrt(n))
# Plot z against exact and reference on the same axes.
# These are analytic densities, not histograms of repeated simulations.

Standardization and normal probability bands

For any variable with finite mean μ and finite positive SD σ, Z = (X−μ)/σ has mean 0 and variance 1. It is standard normal only if X is normal. Using sample summaries gives descriptive standardized values, not automatically a standard-normal reference distribution. For normal X, μ±kσ corresponds to −k≤Z≤k:

Under an absolutely continuous model, interval probability is density area and a point has probability zero. This is not a universal ban on point probabilities: discrete PMFs assign positive masses. Standardization does not eliminate the need to justify a sampling model.

9. Confidence intervals: target, assumptions and coverage

A frequentist confidence interval is a procedure returning a random interval from a sample. Under its model, an exactly calibrated 95% procedure covers the fixed parameter in 95% of repeated samples. An observed interval either contains it or does not; the level is not a posterior probability assigned to μ after seeing the data. Approximate procedures need not attain exact nominal coverage. See NIST’s interpretation.

Known σ and estimated s require different pivots

Here z(p) is the standard-normal p-quantile; t(p, ν) is the p-quantile of Student’s t with ν degrees of freedom. The upper quantile for a two-sided level 1−α is p = 1−α/2.

Independent normal observations, known σ:
Z = (X̄ − μ)/(σ/√n) ~ N(0,1)
CI for μ: x̄ ± z(1−α/2) σ/√n

Independent normal observations, unknown σ, n > 1:
s² = Σ(xᵢ − x̄)²/(n−1)
T = (X̄ − μ)/(S/√n) ~ tₙ₋₁
CI for μ: x̄ ± t(1−α/2, n−1) s/√n

Two-sided z critical values are 1.959964… for 95% and 2.575829… for 99%; the slide’s 1.96 and 2.58 are rounded. With estimated spread under normal sampling, use the t construction, not merely z with s substituted for σ. Outside normal sampling, approximations need justification; neither n = 31 nor a normality pretest guarantees valid inference for dependent days.

import random
from statistics import NormalDist, mean
from math import sqrt

rng = random.Random(20260908)
mu, sigma, n, repeats = 100, 20, 10, 20
samples = [[rng.gauss(mu, sigma) for _ in range(n)] for _ in range(repeats)]
half = NormalDist().inv_cdf(0.975) * sigma/sqrt(n)
intervals = [(mean(x)-half, mean(x)+half) for x in samples]
covered = sum(lo <= mu <= hi for lo, hi in intervals)
print(covered, "of", repeats, "contain the fixed mu")
# The seed is fixed for reproduction, not chosen to force 19 of 20 successes.
# Coverage is a property over repeated samples, not a guarantee for this batch.

Probability bands, mean CIs and predictions have different targets

The prediction interval includes estimation uncertainty and variation of the new observation. Coverage is over joint repetition of the training sample and new observation, not a guarantee that 95% of future values lie in every realized interval. Increasing n does not remove individual variability. A tolerance interval addressing a population proportion is another construction.

import numpy as np
from scipy.stats import norm, t

data = np.array([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], dtype=float)
center, s, n = data.mean(), data.std(ddof=1), len(data)
level = 0.95
zcrit, tcrit = norm.ppf((1+level)/2), t.ppf((1+level)/2, df=n-1)
# Hypothetical known-sigma comparison: set sigma numerically equal to s.
# This does NOT make sigma known for the recorded traffic days.
sigma = s
halves = [zcrit*sigma/np.sqrt(n), tcrit*s/np.sqrt(n), tcrit*s*np.sqrt(1+1/n)]
for target, h in zip(["z mean CI", "t mean CI", "t prediction"], halves):
    print(target, center-h, center+h)
# Formula evaluation only: independent normal sampling is not established.

At the default summaries, the conditional 95% t mean interval is [3993.85, 4522.15]; the hypothetical known-σ z interval is [4004.49, 4511.51]. These calculations are not validated inference about traffic: σ is not known, independent normal sampling is not established, and the 31 days may be the complete descriptive period rather than a random sample from the intended target. Define the target and sampling model first.

Interval–test equivalence requires matching assumptions, statistic, sidedness, level and endpoint convention (section 10). A narrower interval is not automatically more trustworthy: a larger hypothetical n does not repair systematic error, dependence or model misspecification.

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

A null hypothesis H₀ specifies a population/model claim, such as a head probability of 0.5. It is not universally the statement “everything is due to chance.” A test also needs assumptions, a statistic, an alternative, a rule for what counts as at least as extreme, and a significance level chosen before examining the outcome.

A p-value is a probability under the specified null model of outcomes at least as extreme as the observation, using the chosen test’s ordering. It is not the probability that H₀ is true, nor the probability that the result arose “by chance alone.” The NIST testing guide explains the tail-probability interpretation; the ASA’s statement summary emphasizes model assumptions, context and complete reporting.

Exact fair-coin calculation

Let K be the number of heads in ten independent flips with common head probability θ. Under H₀: θ = 0.5, P(K = j) = C(10, j) / 2¹⁰ = C(10, j) / 1024. For eight heads:

• One outcome: P(K = 8) = 45/1024 = 0.0439453125. This is not the test’s p-value.
• Prespecified upper-tail alternative θ > 0.5: P(K ≥ 8) = (45 + 10 + 1)/1024 = 0.0546875.
• Two-sided alternative θ ≠ 0.5: sum outcomes no more probable under H₀ than the observed one. For this symmetric distribution, these are K ≤ 2 or K ≥ 8, giving 112/1024 = 0.109375.

Neither of these p-values rejects at α = 0.05. Comparing just P(K = 8) with α does not implement either test. For a general asymmetric null, do not assume that every two-sided definition is twice the smaller tail. SciPy’s binomtest exposes the alternative explicitly.

A decision is not a proof or an effect size

This chapter uses the rule reject if p ≤ α, otherwise do not reject. Boundary conventions must be consistent; none of the ten-flip p-values equals 0.05. A rejection can be a Type I error. Non-rejection does not prove fairness, equality or equivalence; a small sample may have low power. Statistical significance also does not establish practical importance or a causal explanation.

Confidence sets can be built by inverting a specified family of tests. Then inclusion/exclusion of a null value matches non-rejection/rejection under that same construction, level and boundary convention. An arbitrary interval and a differently defined test need not agree, especially for discrete distributions. A single estimate lying in a central band of raw data is not this test–interval duality.

from scipy.stats import binomtest
for alternative in ['two-sided', 'greater', 'less']:
    result = binomtest(8, n=10, p=0.5, alternative=alternative)
    print(alternative, result.pvalue)
# two-sided 0.109375
# greater   0.0546875
# less      0.9892578125

When comparing forecasting methods, specify the loss, paired observations, time dependence and relevant effect size before selecting a test. A single significance threshold does not establish that two models are equivalent or identify a universally best model.

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

Degrees of freedom depend on the statistic and model, not a universal “sample size minus one” rule. For example, estimating one sample mean leaves n − 1 freely varying residuals; the usual one-sample t statistic has that reference distribution under independent normal observations. The exact binomial test above instead uses n and the null head probability, without a t-style degrees-of-freedom parameter.

Choose a one-sided alternative when a direction is part of the question in advance; use a two-sided alternative for departures in either direction. Selecting the smaller one-sided p-value after observing the result changes the procedure and its error rate.

Power depends on the alternative, sample size, variability and decision rule. It is not determined by α alone, and a composite alternative generally has a power function, not one universal β. Which error is more costly depends on the application; Type I is not inherently the more serious error.

These are exact calculations for independent Bernoulli trials, not a diagnosis of a physical coin or a general forecast-comparison result. Rejection does not logically prove the alternative; non-rejection does not prove equality. Report uncertainty and the experimental design along with the decision.

12. Choosing a test: question, design and assumptions

Start with the quantity to compare or association to describe, the measurement scale, the independent sampling unit, and the null/alternative hypotheses. A list of test names is not an automatic decision rule. In particular, a method that answers a question about means need not be interchangeable with one that answers a question about ranks or distributions.

Parametric does not mean Gaussian; non-parametric does not mean assumption-free

A parametric model specifies a family described by finitely many parameters—for example a normal model or the binomial model in section 10. Non-parametric procedures do not require such a fully specified parametric family, but still have assumptions about sampling, independence, exchangeability, symmetry or the measurement scale, depending on the procedure. Distribution-free null calibration holds under stated conditions, not for arbitrary data. Many, but not all, non-parametric procedures use ranks.

The slides mention non-normal shape, small samples and ordinal measurements as reasons to consider non-parametric methods. These are prompts to examine the design, not sufficient selection rules. A small sample does not remove the assumptions of a rank test, and a normality pretest does not by itself choose the right analysis. There is no universal power ordering: power depends on the alternative, sample size, design and procedure.

Counts and relationships

Paired versus independent comparisons

Pairing comes from meaningful correspondence: two measurements on the same subject, or two model losses on the same held-out case or time point. Equal array lengths alone do not create pairs. For a paired analysis, form differences within the original units; independently sorting the two columns destroys the pairing. Independence between units is a separate condition from dependence within a pair.

The procedures in a row are candidates for different targets and assumptions, not interchangeable buttons. For multiple conditions, plan contrasts and multiplicity handling; an omnibus rejection is not evidence that every pair differs. For ordinal outcomes, use procedures justified for that scale instead of silently treating category codes as equally spaced measurements.

13. Assessing normality: plots, assumptions and tests

Not every parametric method assumes a Gaussian distribution: a binomial model is parametric too. When a method does require normality, identify which quantity the assumption concerns—for example model errors or paired differences, not automatically every raw variable. Graphs and tests provide complementary diagnostics; neither certifies a model from a finite sample.

1. Summary checks are not a calibrated normality test

The slides propose comparing mean and median within 1% of the range and checking the 68–95–99.7 percentages with at least 50 observations. Treat those as informal teaching heuristics, not a universal threshold or required minimum. Symmetry alone does not imply normality (a uniform distribution is symmetric too). A random normal sample need not meet exact population percentages or have equal sample mean and median. Outliers can also inflate the range and make a range-based tolerance misleading.

2. Histograms reveal shape, but depend on the bins

Look for asymmetry, multiple modes and unusual tails, checking more than one sensible binning. Small samples make shape particularly uncertain. If a normal density curve is overlaid, use a density-normalized histogram or scale the curve to expected bin counts; a probability density and an unscaled count histogram have different vertical units. See the histogram widget in section 2 for how counts change with bin boundaries.

3. A Q–Q plot pairs ordered observations with theoretical quantiles

Sort the observations, assign increasing plotting probabilities pᵢ, then calculate zᵢ = Φ⁻¹(pᵢ), where Φ is the standard normal cumulative distribution function. Plot zᵢ horizontally and the ordered observations vertically. A roughly linear pattern is a useful diagnostic; a normal sample does not give an exact line. The NIST normal probability plot guide uses a different plotting-position convention (approximate order-statistic medians). Here we retain the slides’ midpoint convention, pᵢ = (i − ½)/n. These are chosen probability positions, not a claim that actual random observations enclose exactly equal probability mass.

The seven original values are −4, −3, 0.8, 1.8, 3.9, 6.2 and 6.5. Their mean is 1.742857. Define dₙ = √[Σ(xᵢ − x̄)²/n] = 3.845538: this reproduces the slides’ rounded 3.85. The n−1 sample standard deviation is instead s = 4.153656. This drawing deliberately uses dₙ, so the reference line is y = x̄ + dₙz; it is a descriptive fitted line, not an uncertainty interval.

The final column is qᵢ = x̄ + dₙzᵢ in the observations’ units. If qᵢ were used on the horizontal axis instead, the reference would be y = x, visually 45° only with equal unit scaling. Do not confuse these two axis conventions. Table and figure use unrounded pᵢ throughout; the displayed 0.071429 is not substituted back into the inverse CDF.

from statistics import mean, pstdev, NormalDist
x = sorted([-4, -3, 0.8, 1.8, 3.9, 6.2, 6.5])
mu, d = mean(x), pstdev(x)  # denominator n, matching this figure
for i, observed in enumerate(x, 1):
    p = (i - 0.5) / len(x)
    z = NormalDist().inv_cdf(p)
    print(i, observed, p, z, mu + d * z)

Software defaults are part of the specification. In statsmodels.qqplot, the default offset is a = 0, producing i/(n+1), and line='q' draws a quartile-based line. Those defaults do not reproduce this table or its mean/scale line. Use a = 0.5 for the midpoint positions and explicitly supply the intended reference line when recreating this figure.

4. Statistical tests answer a specified null hypothesis

Shapiro–Wilk tests the normal family with unspecified location and scale; it is not simply a KS test with fitted parameters. Anderson–Darling and D’Agostino’s K-squared are other normality procedures with their own calibration and requirements. A chi-square goodness-of-fit test requires choices of bins and appropriate expected counts. There is no universal power ranking: small samples can miss meaningful departures, while large samples can detect departures too small to matter for the intended analysis.

The usual one-sample Kolmogorov–Smirnov test instead requires a fully specified continuous null distribution. Estimating its parameters on the tested sample invalidates the usual null calibration; a suitable fitted-parameter procedure such as the normal Lilliefors test, or a correctly calibrated simulation procedure, is needed. Independence and the sampling/model assumptions still matter.

"""Same seeded sample, three explicitly different null hypotheses."""
import numpy as np
from scipy.stats import shapiro, kstest

# Preserve the lecture's legacy random stream, not default_rng's different stream.
data = 20 * np.random.RandomState(20).randn(100) + 100
alpha = 0.05
tests = [
    ("Some normal distribution", shapiro(data)),
    ("Specified N(0, 1)", kstest(data, "norm")),
    ("Specified N(100, 20²)", kstest(data, "norm", args=(100, 20))),
]
for null, result in tests:
    decision = "Reject H0" if result.pvalue <= alpha else "Do not reject H0"
    print(f"{null}: statistic={result.statistic:.6g}, p={result.pvalue:.6g}; {decision}")
# The last null uses the known generating parameters, not estimates from data.
# Non-rejection is not proof of normality; p printed as 0 can reflect numerical limits.

This is synthetic data generated with mean 100 and standard deviation 20, not an observed data set known to be normal. A rejection of the standard N(0,1) null is not evidence against every normal distribution. The third test uses the generating parameters fixed independently of the sample. A p-value printed as zero can reflect numerical limits, not a literal probability of zero. Report the null, assumptions, statistic and p-value; say “do not reject” rather than “accept” or “proved Gaussian.”

14. Student’s t, Diebold–Mariano and forecast alignment

Student’s t is a family indexed by degrees of freedom. For an independent normal sample with unknown variance, the standardized sample mean uses a t reference with n−1 degrees of freedom. Unknown variance alone does not make every statistic exactly t-distributed. Other uses, including the modified forecast-comparison test here, have their own justification and limitations.

Define the loss and the comparison before testing

For aligned forecasts of the same target, let dₜ = L(e₁ₜ) − L(e₂ₜ). The Diebold–Mariano null is equal expected loss, E[dₜ] = 0, for the chosen loss and forecast setting—not equal predictions at every time or equal accuracy under every metric. Positive average d means the first forecast has larger observed loss. The author’s retrospective working paper emphasizes comparing forecasts and making assumptions about their loss differences, rather than requiring knowledge of a model’s architecture.

The unmodified statistic is d̄ / √(Ω̂/n), where Ω̂ estimates the long-run variance of d. Its normal reference is asymptotic and needs suitable stationarity, moment and dependence conditions plus consistent variance estimation. Serial dependence, nonstationarity, parameter estimation and nested-model settings require care. “Model-free” does not mean assumption-free.

The printed statistic uses a modified test, not a normal critical value

The implementation linked by the slides applies the Harvey–Leybourne–Newbold adjustment and evaluates a two-sided t tail with n−1 degrees of freedom. The forecast package also documents a modified test. For this example n = 12 and h = 1: γ̂₀ = Σ(dₜ − d̄)²/n, the unmodified statistic is d̄/√(γ̂₀/n), and the adjustment multiplies it by √(11/12). This equals the ordinary t statistic of the loss differences algebraically; that identity does not prove those differences independent or normal.

Use the matching 5% two-sided critical values ±2.200985 from t₁₁, not ±1.96. With h = 1 this implementation keeps only lag zero; the horizon does not establish that all other autocovariances vanish. The modified t reference is not an exact finite-sample guarantee for arbitrary forecast errors, especially with only twelve observations.

Audit the two conflicting MLP sequences

The slide 66 list begins 414.231, 417.662, 381.815; the slide 67 executable example moves 414.231 to the end. Actual and LSTM arrays are unchanged. Without timestamped prediction exports and training/evaluation provenance, these sources do not identify which alignment is correct. Do not rotate predictions to obtain a more favorable score. Both calculations below preserve their respective printed orderings.

The code ordering reproduces modified DM = 2.933300 and p = 0.01360984. The list ordering instead gives 3.429053 and 0.00563247. Both arithmetic decisions reject at 5% under this specified procedure, but neither establishes that the data preparation or forecast alignment is valid. The p-value is not the probability that the difference was “due to chance,” and rejection is not proof of universal LSTM superiority.

For comparison, the unmodified code-order statistic is 3.063732; a normal-tail calculation for that statistic gives 0.00218594. Do not combine a statistic from one variant with a p-value or critical value from another. The LSTM values differ from all twelve actual values by at most 0.001; their mean squared error is 5.000000e-7. This calls for checking the prediction pipeline and data availability, not assuming either extraordinary generalization or leakage from these numbers alone.

"""Reproduce the two printed orderings; this does not validate forecast provenance."""
import json
import math
from scipy.stats import t, norm

ACTUAL = [417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432]
MLP_LIST = [414.231, 417.662, 381.815, 423.952, 469.055, 471.502,
            547.532, 639.609, 596.040, 477.471, 445.010, 368.674]
MLP_CODE = [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]


def calculate(actual, first, second):
    # Deliberately limited to this 12-position, squared-loss, h=1 worked example.
    if any(not isinstance(x, list) or len(x) != 12 for x in (actual, first, second)):
        raise ValueError("This example requires three aligned lists of length 12")
    if any(type(v) not in (int, float) or not math.isfinite(v) or abs(v) > 1e6
           for x in (actual, first, second) for v in x):
        raise ValueError("Expected finite numeric values of magnitude at most 1e6")
    n = len(actual)
    loss1 = [(a - p)**2 for a, p in zip(actual, first)]
    loss2 = [(a - p)**2 for a, p in zip(actual, second)]
    differences = [a - b for a, b in zip(loss1, loss2)]
    mean = math.fsum(differences) / n
    gamma0 = math.fsum((d - mean)**2 for d in differences) / n
    if gamma0 <= 0:
        raise ValueError("A positive loss-difference variance is required")
    raw = mean / math.sqrt(gamma0 / n)
    factor = math.sqrt((n - 1) / n)  # Harvey–Leybourne–Newbold factor at h=1
    modified = factor * raw
    return dict(mean=mean, mse1=math.fsum(loss1)/n, mse2=math.fsum(loss2)/n,
                differences=differences, gamma0=gamma0, raw=raw, factor=factor,
                statistic=modified, pvalue=float(2*t.sf(abs(modified), n-1)),
                critical=float(t.ppf(0.975, n-1)), normalP=float(2*norm.sf(abs(raw))))


def example():
    return dict(n=12, h=1, loss="squared error", df=11, alpha=0.05,
                actual=ACTUAL, lstm=LSTM, mlpList=MLP_LIST, mlpCode=MLP_CODE,
                listResult=calculate(ACTUAL, MLP_LIST, LSTM),
                codeResult=calculate(ACTUAL, MLP_CODE, LSTM))


if __name__ == "__main__":
    print(json.dumps(example(), indent=2, allow_nan=False))

Critical-difference diagrams summarize a defined analysis

A critical-difference diagram can show average ranks alongside pairwise post-hoc results. A crossbar denotes non-rejection under the chosen procedure and threshold, not equivalence or identical performance. Non-rejection need not be transitive, and ranks alone do not determine adjusted p-values. Distinguish a common critical-distance procedure from pairwise Wilcoxon tests with Holm adjustment; do not manufacture one shared distance for arbitrary pairwise results.

The slides describe five labelled methods, but the available extraction does not provide the full score matrix or adjusted pairwise p-values. Its specific crossbars cannot be numerically reproduced from this material alone. For multiple forecasting methods, choose valid comparisons and multiplicity handling for the actual dependence structure; adding models does not automatically justify a repeated-measures test on overlapping fold scores.

Check your understanding