- Compute the mean, median, standard deviation and IQR, and detect outliers
- Find probabilities and quantiles of the normal distribution with
scipy.stats.norm - Build a confidence interval for a mean and interpret a t-test's p-value correctly
- Compute Pearson and Spearman correlation and tell correlation from causation
A food delivery service in Baku promises delivery “in 15 minutes on average”. Murad timed 10 orders. Is that enough to check the promise? If one order took 45 minutes, how does it affect the result? Statistics answers in three steps: describe the data, model it with a distribution and infer from the sample to the whole population. In Python, NumPy and scipy.stats are all you need.
Descriptive statistics and outliers
Measures of centre give a “typical” value: the mean (arithmetic average), the median (the middle of the sorted data) and the mode (the most frequent value). Measures of spread show how scattered the values are: the variance, the standard deviation and the interquartile range IQR = Q3 − Q1. For a sample we divide the variance by n − 1 (Bessel's correction), which makes the estimate unbiased.
- nsample size
- s²sample variance (its unit is the square of the data's unit)
- ssample standard deviation (same unit as the data); in NumPy
std(ddof=1)
import numpy as np
from scipy import stats
minutes = np.array([12, 15, 11, 14, 13, 45, 13, 16, 14, 13])
print('mean :', minutes.mean())
print('median:', np.median(minutes))
print('mode :', stats.mode(minutes, keepdims=False).mode)
print('var :', round(minutes.var(ddof=1), 2))
print('std :', round(minutes.std(ddof=1), 2))
q1, q3 = np.percentile(minutes, [25, 75])
iqr = q3 - q1
print('Q1, Q3, IQR:', q1, q3, iqr)
low, high = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print('outliers:', minutes[(minutes < low) | (minutes > high)])▸ Expected output
mean : 16.6 median: 13.5 mode : 13 var : 101.6 std : 10.08 Q1, Q3, IQR: 13.0 14.75 1.75 outliers: [45]
- Q1, Q3first and third quartile (25% and 75%)
- IQRQ3 − Q1, the width of the middle half of the data
Tukey's rule: values outside this interval count as outliers. Here: 14.75 + 1.5 · 1.75 = 17.375, so 45 is an outlier.
The normal distribution
Quantities that are the sum of many small independent effects (measurement errors, test scores, heights) are often close to a normal distribution: a bell-shaped curve, symmetric around the mean. It is defined by two parameters: the mean μ and the standard deviation σ. Any value can be standardised with a z-score, which shows how many σ it lies from the mean.
- μmean of the distribution
- σstandard deviation (σ > 0)
- f(x)probability density; P(a < X < b) is the area under the curve between a and b
- zz-score: the position in the standard normal distribution (μ = 0, σ = 1)
scipy.stats.norm creates a distribution object: cdf(x) = P(X ≤ x) (the cumulative distribution function), sf(x) = P(X > x) = 1 − cdf, and ppf(q) is the inverse — the quantile for a given probability. If test scores are normal with μ = 70 and σ = 10:
from scipy import stats
scores = stats.norm(loc=70, scale=10)
print(round(scores.cdf(85), 4))
print(round(scores.sf(85), 4))
print(round(scores.cdf(80) - scores.cdf(60), 4))
print(round(scores.ppf(0.90), 2))
z = (85 - 70) / 10
print(z, round(stats.norm.cdf(z), 4))▸ Expected output
0.9332 0.0668 0.6827 82.82 1.5 0.9332
Test scores follow N(70, 10²). a) What is the probability that a random student scores more than 85? b) What score puts you in the top 10%? Φ(1.5) = 0.9332, z₀.₉₀ = 1.2816.
Show solutionHide solution
P(X > 85) = 1 − Φ(1.5) = 1 − 0.9332 = 0.0668, about 6.7%.
b) x = μ + z₀.₉₀ · σ = 70 + 1.2816 · 10 ≈ 82.8 points.
Both answers match the code (
sf(85) and ppf(0.90)).Samples, the central limit theorem and the standard error
We usually see not the whole population (all orders) but only a sample. Each new sample has a slightly different mean. The central limit theorem says: when n is large enough, the distribution of sample means becomes close to normal even if the original distribution is skewed, and its standard deviation — the standard error — equals σ/√n. Below we take 5000 samples of size 30 from a skewed exponential distribution:
- SEstandard error of the mean
- σ, spopulation and sample standard deviation
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(1)
population = rng.exponential(scale=10, size=100_000)
means = rng.choice(population, size=(5000, 30)).mean(axis=1)
print(f'population: mean {population.mean():.2f}, std {population.std():.2f}')
print(f'sample means: mean {means.mean():.2f}, std {means.std():.2f}')
print(f'sigma / sqrt(n) = {population.std() / np.sqrt(30):.2f}')
fig, (a, b) = plt.subplots(1, 2, figsize=(8, 3), layout='constrained')
a.hist(population, bins=50)
a.set_title('Population (skewed)')
b.hist(means, bins=40)
b.set_title('Means of samples, n = 30')
plt.show()▸ Expected output
population: mean 9.96, std 9.92 sample means: mean 9.94, std 1.81 sigma / sqrt(n) = 1.81
Confidence intervals and the t-test
When σ is unknown we replace it with s and use Student's t-distribution with n − 1 degrees of freedom instead of the normal one — for small samples its “tails” are heavier. Ten 100-gram packs of tea from a factory line were weighed:
- t*critical value of the t-distribution; α = 0.05 for 95%
- n − 1degrees of freedom
A 95% confidence interval for the mean. stats.t.ppf(0.975, df=n - 1) returns t*, and stats.t.interval returns the whole interval.
import numpy as np
from scipy import stats
packs = np.array([99.2, 100.4, 98.7, 100.9, 99.5, 100.1, 99.8, 98.9, 100.6, 99.4])
n = packs.size
mean = packs.mean()
s = packs.std(ddof=1)
se = s / np.sqrt(n)
t_crit = stats.t.ppf(0.975, df=n - 1)
print(f'mean = {mean:.2f} g, s = {s:.3f} g, SE = {se:.3f} g')
print(f't* = {t_crit:.3f}')
print(f'95% CI: [{mean - t_crit * se:.2f}; {mean + t_crit * se:.2f}] g')
low, high = stats.t.interval(0.95, df=n - 1, loc=mean, scale=se)
print(round(low, 2), round(high, 2))▸ Expected output
mean = 99.75 g, s = 0.738 g, SE = 0.233 g t* = 2.262 95% CI: [99.22; 100.28] g 99.22 100.28
The exam preparation time of 16 students: x̄ = 52 hours, s = 8 hours. Find a 95% confidence interval for the mean. t₀.₉₇₅(15) = 2.131.
Show solutionHide solution
Margin of error = t* · SE = 2.131 · 2 = 4.262.
Interval: 52 ± 4.26 → [47.74, 56.26] hours.
With 64 students the SE would drop to 1 and the interval would become about twice as narrow (the √n law).
A hypothesis test asks: if the null hypothesis H₀ (“there is no difference”) were true, how likely would a difference as large as ours or larger be? That probability is the p-value. If p < α (usually 0.05), we reject H₀. A one-sample t-test compares a mean with a given number, and a two-sample test compares the means of two groups; since the groups' variances may differ, we choose Welch's version (equal_var=False).
- x̄₁, x̄₂group means
- s₁², s₂²sample variances of the groups
- n₁, n₂group sizes
Welch's t statistic: the difference of the means divided by its standard error. The larger |t|, the smaller p.
import numpy as np
from scipy import stats
packs = np.array([99.2, 100.4, 98.7, 100.9, 99.5, 100.1, 99.8, 98.9, 100.6, 99.4])
res1 = stats.ttest_1samp(packs, popmean=100)
print(f'one-sample: t = {res1.statistic:.3f}, p = {res1.pvalue:.3f}')
group_a = np.array([72, 85, 78, 90, 66, 81, 77, 88, 74, 83])
group_b = np.array([80, 91, 86, 94, 79, 88, 85, 97, 82, 90])
print(group_a.mean(), group_b.mean())
res2 = stats.ttest_ind(group_b, group_a, equal_var=False)
print(f'Welch: t = {res2.statistic:.3f}, p = {res2.pvalue:.4f}')▸ Expected output
one-sample: t = -1.071, p = 0.312 79.4 87.2 Welch: t = 2.581, p = 0.0194
Correlation
The Pearson correlation coefficient r measures the strength and direction of a linear relationship between two variables, from −1 to 1; r² shows what share of the variation in y is explained linearly by x. The Spearman coefficient ρ does the same calculation on ranks, so it is robust to outliers and captures any monotonic relationship.
- r−1 ≤ r ≤ 1; the sign gives the direction, |r| the strength
- x̄, ȳthe mean of each variable
import numpy as np
from scipy import stats
hours = np.array([1, 2, 2, 3, 4, 5, 5, 6, 7, 8])
score = np.array([52, 55, 61, 60, 68, 70, 75, 74, 82, 88])
r, p = stats.pearsonr(hours, score)
rho, _ = stats.spearmanr(hours, score)
print(f'Pearson r = {r:.3f}, p = {p:.1e}')
print(f'Spearman rho = {rho:.3f}')
print(f'r squared = {r ** 2:.3f}')
print(np.corrcoef(hours, score).round(3))▸ Expected output
Pearson r = 0.980, p = 6.5e-07 Spearman rho = 0.957 r squared = 0.961 [[1. 0.98] [0.98 1. ]]
np.corrcoef returns the correlation matrix.x = (1, 2, 3, 4, 5), y = (2, 4, 5, 4, 5). Find the Pearson correlation coefficient.
Show solutionHide solution
xᵢ − x̄: −2, −1, 0, 1, 2; yᵢ − ȳ: −2, 0, 1, 0, 1.
Sum of products: 4 + 0 + 0 + 0 + 2 = 6.
∑ (xᵢ − x̄)² = 10, ∑ (yᵢ − ȳ)² = 4 + 0 + 1 + 0 + 1 = 6.
r = 6 / √(10 · 6) = 6 / √60 ≈ 0.775; r² ≈ 0.6.
Check:
np.corrcoef(x, y)[0, 1] → 0.7746.For the data array, print on the first line the mean (1 decimal place), the median and the sample standard deviation (ddof=1, 2 decimal places), separated by spaces. On the second line print the list of outliers by Tukey's 1.5 · IQR rule.
import numpy as np
data = np.array([48, 52, 50, 47, 53, 51, 49, 95, 50, 52])
# line 1: mean (1 decimal), median, sample std (2 decimals)
# line 2: outliers as a list▸ Expected output
54.7 50.5 14.28 [95]
A lab made 8 measurements. Print the sample mean on the first line and, on the second line, the lower and upper bounds of the 95% confidence interval for the mean (separated by a space); round everything to 2 decimal places. Find t* with stats.t.ppf.
import numpy as np
from scipy import stats
sample = np.array([23.1, 24.5, 22.8, 25.0, 23.9, 24.2, 23.5, 24.8])
# mean, then the 95% confidence interval▸ Expected output
23.98 23.31 24.64
Key points
- Sample variance divides by n − 1 (
ddof=1); for skewed data the median and IQR are more reliable. - Tukey's rule: values outside [Q1 − 1.5 · IQR, Q3 + 1.5 · IQR] are outliers.
- z = (x − μ) / σ;
norm.cdf,norm.sfandnorm.ppfgive probabilities and quantiles; the 68–95–99.7 rule. - Standard error SE = s / √n; the interval for a mean is x̄ ± t* · SE.
- Reject H₀ when p < α, but p is neither the probability of H₀ nor the size of the effect.
- Pearson's r measures linear relationships and Spearman's ρ works with ranks; correlation does not prove causation.
Check yourself
10 questions. Every correct answer earns XP.