Statistics for Python Intermediate¶
When you'd use this
scipy.stats, distributions, hypothesis testing, correlation and regression.
Summarize and reason about data: distributions, hypothesis tests, correlation — for analysis and validating experiments.
Descriptive statistics¶
Summarize data with mean, median, variance, and quantiles — the first look at any dataset.
import numpy as np
from scipy import stats
data = np.array([23, 45, 12, 67, 34, 89, 56, 78, 43, 21, 55, 67, 34, 90, 12])
# Central tendency
print(f"Mean: {np.mean(data):.2f}") # 48.40
print(f"Median: {np.median(data):.2f}") # 45.00
print(f"Mode: {stats.mode(data, keepdims=True).mode[0]}") # 12 or 34 or 67
# Spread
print(f"Std: {np.std(data, ddof=1):.2f}") # 25.88 (sample std)
print(f"Var: {np.var(data, ddof=1):.2f}") # 669.83
print(f"Range: {np.ptp(data)}") # 78 (max - min)
print(f"IQR: {stats.iqr(data):.2f}") # 35.50
# Percentiles
print(f"25th: {np.percentile(data, 25):.2f}") # 23.00
print(f"75th: {np.percentile(data, 75):.2f}") # 67.00
print(f"90th: {np.percentile(data, 90):.2f}") # 84.80
# Shape
print(f"Skewness: {stats.skew(data):.4f}") # ~0 for symmetric
print(f"Kurtosis: {stats.kurtosis(data):.4f}") # ~0 for normal-like
Probability distributions¶
Model randomness with common distributions for sampling and inference.
from scipy.stats import norm, t, chi2, binom, poisson, expon
# ─── Normal distribution ──────────────────────────
# PDF — probability density at a point
print(norm.pdf(0, loc=0, scale=1)) # 0.3989 (peak of standard normal)
# CDF — probability of X ≤ x
print(norm.cdf(1.96)) # 0.975 (97.5% below 1.96σ)
print(norm.cdf(0)) # 0.5 (50% below mean)
# Inverse CDF (quantile function)
print(norm.ppf(0.975)) # 1.96
print(norm.ppf(0.95)) # 1.645
# Random samples
samples = norm.rvs(loc=100, scale=15, size=1000, random_state=42)
print(f"Sample mean: {samples.mean():.2f}, std: {samples.std():.2f}")
# ─── t-distribution ──────────────────────────────
# Used when sample size is small and population std is unknown
t_critical = t.ppf(0.975, df=29) # 95% CI, 30 samples
print(f"t-critical (df=29): {t_critical:.4f}") # 2.045
# ─── Binomial ────────────────────────────────────
# P(X = k) for n trials with probability p
print(binom.pmf(k=7, n=10, p=0.5)) # P(7 heads in 10 flips)
print(binom.cdf(k=7, n=10, p=0.5)) # P(≤ 7 heads)
# ─── Poisson ─────────────────────────────────────
# Events per unit time (λ = average rate)
print(poisson.pmf(k=5, mu=3)) # P(5 events when average is 3)
Hypothesis testing¶
Decide whether an observed effect is real or chance with t-tests and friends.
One-sample t-test¶
Test whether a sample's mean differs from a known reference value — e.g. "is average height different from 170cm?"
# Question: Is the average height different from 170cm?
heights = np.array([172, 168, 175, 171, 169, 174, 167, 173, 170, 176,
171, 168, 172, 174, 169, 173, 175, 170, 168, 172])
t_stat, p_value = stats.ttest_1samp(heights, popmean=170)
print(f"t-statistic: {t_stat:.4f}") # 2.456
print(f"p-value: {p_value:.4f}") # 0.0237
alpha = 0.05
if p_value < alpha:
print("Reject H0: Average height IS different from 170cm")
else:
print("Fail to reject H0: No significant difference")
Two-sample t-test¶
Compare the means of two independent groups (or paired before/after measurements) to see if they differ significantly.
# Question: Do two groups have different means?
group_a = np.array([85, 90, 78, 92, 88, 76, 95, 89, 91, 84])
group_b = np.array([75, 80, 72, 85, 78, 70, 82, 77, 79, 74])
t_stat, p_value = stats.ttest_ind(group_a, group_b)
print(f"t-statistic: {t_stat:.4f}") # 3.826
print(f"p-value: {p_value:.6f}") # 0.001229
# For paired samples (before/after):
t_stat, p_value = stats.ttest_rel(group_a, group_b)
Chi-squared test (independence)¶
Test whether two categorical variables are related — e.g. "does product preference depend on gender?" — from a contingency table of counts.
# Question: Is there a relationship between gender and product preference?
observed = np.array([
[45, 30, 25], # Male: Product A, B, C
[35, 40, 25], # Female: Product A, B, C
])
chi2_stat, p_value, dof, expected = stats.chi2_contingency(observed)
print(f"χ² statistic: {chi2_stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Degrees of freedom: {dof}")
print(f"Expected frequencies:\n{expected}")
ANOVA (multiple groups)¶
Compare the means of three or more groups at once — avoids the inflated error rate of running many pairwise t-tests.
# Question: Do 3+ groups have different means?
group1 = [85, 90, 78, 92, 88]
group2 = [75, 80, 72, 85, 78]
group3 = [92, 95, 88, 97, 91]
f_stat, p_value = stats.f_oneway(group1, group2, group3)
print(f"F-statistic: {f_stat:.4f}")
print(f"p-value: {p_value:.6f}")
Correlation¶
Measure how two variables move together — but remember correlation isn't causation.
import numpy as np
from scipy import stats
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.array([2.1, 3.8, 6.2, 7.9, 10.1, 12.3, 13.8, 16.1, 18.0, 19.9])
# Pearson (linear correlation)
r, p_value = stats.pearsonr(x, y)
print(f"Pearson r: {r:.4f}, p-value: {p_value:.6f}")
# r ≈ 0.9997 (nearly perfect linear relationship)
# Spearman (rank correlation — handles non-linear monotonic)
rho, p_value = stats.spearmanr(x, y)
print(f"Spearman ρ: {rho:.4f}")
# Correlation matrix (with pandas)
import pandas as pd
df = pd.DataFrame({"x": x, "y": y, "z": x**2})
print(df.corr())
# x y z
# x 1.000000 0.999720 0.974679
# y 0.999720 1.000000 0.970818
# z 0.974679 0.970818 1.000000
Linear regression¶
Fit a line to model and predict a numeric outcome from features.
from scipy.stats import linregress
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.array([2.5, 4.1, 5.8, 8.2, 9.5, 12.1, 13.7, 15.9, 18.2, 19.8])
result = linregress(x, y)
print(f"Slope: {result.slope:.4f}") # ~1.95
print(f"Intercept: {result.intercept:.4f}") # ~0.57
print(f"R-squared: {result.rvalue**2:.4f}") # ~0.998
print(f"p-value: {result.pvalue:.2e}") # very small
print(f"Std error: {result.stderr:.4f}")
# Predict
x_new = 12
y_pred = result.slope * x_new + result.intercept
print(f"Prediction for x={x_new}: {y_pred:.2f}")
# Plot
import matplotlib.pyplot as plt
plt.scatter(x, y, label="Data")
plt.plot(x, result.slope * x + result.intercept, "r-", label=f"Fit (R²={result.rvalue**2:.3f})")
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("Linear Regression")
plt.show()
Confidence intervals¶
Quantify uncertainty around an estimate instead of a single point.
# CI for the mean
data = np.array([23, 25, 28, 30, 26, 24, 27, 29, 31, 22])
n = len(data)
mean = np.mean(data)
se = stats.sem(data) # standard error of the mean
# 95% confidence interval
ci = stats.t.interval(confidence=0.95, df=n-1, loc=mean, scale=se)
print(f"Mean: {mean:.2f}")
print(f"95% CI: ({ci[0]:.2f}, {ci[1]:.2f})")
# Interpretation: We're 95% confident the true mean is in this interval
# Bootstrap confidence interval (non-parametric)
def bootstrap_ci(data, n_bootstrap=10000, confidence=0.95):
rng = np.random.default_rng(42)
means = [rng.choice(data, size=len(data), replace=True).mean()
for _ in range(n_bootstrap)]
lower = np.percentile(means, (1 - confidence) / 2 * 100)
upper = np.percentile(means, (1 + confidence) / 2 * 100)
return lower, upper
ci_boot = bootstrap_ci(data)
print(f"Bootstrap 95% CI: ({ci_boot[0]:.2f}, {ci_boot[1]:.2f})")
A/B Testing¶
Compare two variants statistically to decide if a change actually helped.
# Conversion rates: variant A vs variant B
visitors_a, conversions_a = 10000, 350 # 3.5%
visitors_b, conversions_b = 10000, 400 # 4.0%
# Two-proportion z-test
from statsmodels.stats.proportion import proportions_ztest
count = np.array([conversions_a, conversions_b])
nobs = np.array([visitors_a, visitors_b])
z_stat, p_value = proportions_ztest(count, nobs, alternative="two-sided")
print(f"z-statistic: {z_stat:.4f}")
print(f"p-value: {p_value:.4f}")
if p_value < 0.05:
print("Statistically significant difference!")
print(f"Variant B converts {(400/10000 - 350/10000)*100:.1f}% more")
Practice Exercises¶
- Generate 1000 samples from a normal distribution and verify the 68-95-99.7 rule.
- Perform a t-test on two groups of exam scores and interpret the result.
- Compute correlation between all pairs of features in a dataset and create a heatmap.
- Implement bootstrap for estimating the confidence interval of the median.
- Design and analyze an A/B test — compute sample size needed for 80% power.
- Fit a linear regression and plot residuals to check assumptions.
💬 Discussion
Have a question about this topic? Found an error? Share your thoughts below.