Statistics for Python: A Practical Guide to Statistical Analysis

The Python Statistics Toolkit

Python does not have a single statistics library. Instead, you draw from several packages, each with a different purpose. Knowing which to reach for first saves a lot of time.

Quick-Reference Function Table

Setting Up Python for Statistics

You can run Python statistics in a local environment or in the browser. Google Colab requires no installation and gives you a Jupyter-style notebook immediately. For local work, install Python 3.9 or later and the libraries below.

pip install numpy pandas scipy statsmodels scikit-learn matplotlib seaborn

For reproducible projects, create a virtual environment first so library versions stay consistent:

python -m venv stats-env
source stats-env/bin/activate          # macOS / Linux
stats-env\Scripts\activate             # Windows
pip install numpy pandas scipy statsmodels matplotlib seaborn

Then start a Jupyter notebook:

pip install jupyterlab
jupyter lab

Standard Imports

Most statistical analysis notebooks begin with these imports. The aliases (np, pd, etc.) are conventional and widely understood.

import numpy as np                   # numerical arrays
import pandas as pd                  # dataframes and series
import scipy.stats as stats           # distributions and tests
import statsmodels.api as sm          # regression and modeling
import statsmodels.formula.api as smf  # formula-style regression
import matplotlib.pyplot as plt       # plotting
import seaborn as sns                 # statistical visualization

The Working Dataset

All examples in this guide use one dataset: student exam scores alongside study hours and course group. This keeps the numbers consistent across descriptive statistics, correlation, hypothesis testing, and regression.

data = pd.DataFrame({
    "hours_studied": [2, 4, 5, 6, 7, 5, 7, 8, 9, 10],
    "exam_score":   [55, 62, 68, 71, 75, 72, 80, 85, 88, 93],
    "group":         ["Control"]*5 + ["Tutored"]*5
})

print(data)

Descriptive Statistics in Python

Descriptive statistics summarize what is in your data. They do not draw conclusions about populations or test hypotheses — they describe the sample you have. The main measurements cover center, spread, and shape.

Mean, Median, and Mode

The mean is the arithmetic average. The median is the middle value when data are sorted. The mode is the value that appears most often.

scores = data["exam_score"]

mean   = scores.mean()      # 74.9
median = scores.median()    # 73.5
mode   = scores.mode()      # pandas returns a Series

print(f"Mean:   {mean}")
print(f"Median: {median}")
print(f"Mode:   {mode.values}")

Mean:   74.9
Median: 73.5
Mode:   [55 62 68 71 72 75 80 85 88 93]  (all values appear once)

The mean and median differ when data are skewed. Here they are close (74.9 vs 73.5), which is consistent with the roughly symmetric spread of the scores. When a dataset has extreme high values, the mean rises above the median — income data is a classic example. The mode tells you the most common value, but in continuous measurement it is often not very informative.

Standard Deviation and Variance — The Critical ddof Distinction

This is the single most common Python statistics mistake. NumPy and pandas calculate standard deviation differently by default.

The sample variance formula divides by n−1 to correct for the fact that sample means are closer to sample values than the true population mean is. This correction — called Bessel's correction — makes the sample variance an unbiased estimator of the population variance.

scores = data["exam_score"]

# Sample statistics (ddof=1) — use for most analyses
sample_std = scores.std()            # pandas default: ddof=1
sample_var = scores.var()            # pandas default: ddof=1

# Population statistics (ddof=0)
pop_std = np.std(scores)             # numpy default: ddof=0
pop_var = np.var(scores)             # numpy default: ddof=0

# Explicit is clearest
sample_std_explicit = np.std(scores, ddof=1)   # same as pandas
pop_std_explicit    = np.std(scores, ddof=0)   # same as default numpy

data_range = scores.max() - scores.min()     # 93 - 55 = 38

print(f"Sample SD (ddof=1):     {sample_std:.2f}")    # 13.21
print(f"Population SD (ddof=0): {pop_std:.2f}")       # 12.54
print(f"Range: {data_range}")

Sample SD (ddof=1):     13.21
Population SD (ddof=0): 12.54
Range: 38

Quartiles and IQR

The interquartile range (IQR) is Q3 minus Q1. It measures the middle 50% of the data and is more resistant to extreme values than the standard deviation.

q1  = scores.quantile(0.25)    # 67.25
q2  = scores.quantile(0.50)    # 73.5
q3  = scores.quantile(0.75)    # 83.25
iqr = q3 - q1                    # 16.0

print(f"Q1: {q1}, Median: {q2}, Q3: {q3}, IQR: {iqr}")

# Outlier detection using IQR rule
lower_fence = q1 - 1.5 * iqr    # 43.25
upper_fence = q3 + 1.5 * iqr    # 107.25
outliers = scores[(scores < lower_fence) | (scores > upper_fence)]
print(f"Outliers: {outliers.values}")  # none in this dataset

The describe() Function

df.describe() returns a convenient summary but is not a complete statistical analysis. It gives you count, mean, standard deviation (using ddof=1), min, quartiles, and max for every numeric column.

print(data.describe())

       hours_studied  exam_score
count      10.000000   10.000000
mean        6.300000   74.900000
std         2.406011   13.208524
min         2.000000   55.000000
25%         5.000000   67.250000
50%         6.500000   73.500000
75%         8.250000   83.250000
max        10.000000   93.000000

The std row here uses ddof=1, consistent with pandas defaults. Think of describe() as a starting point, not an endpoint.

Probability Distributions in Python

A probability distribution describes how likely different values are. SciPy's stats module covers dozens of distributions, each with consistent methods: pdf() for the probability density function, cdf() for the cumulative distribution function, ppf() for the inverse CDF (also called the percent point function), and rvs() to draw random samples.

Normal Distribution

The normal distribution is bell-shaped and fully described by its mean μ and standard deviation σ. A z-score converts any value to units of standard deviations from the mean.

from scipy import stats

mu, sigma = 74.9, 13.21

# Probability density at a point (not a probability for continuous data)
density_at_75 = stats.norm.pdf(75, loc=mu, scale=sigma)

# Cumulative probability: P(X <= 70)
prob_below_70 = stats.norm.cdf(70, loc=mu, scale=sigma)

# Inverse CDF: what score marks the 90th percentile?
score_90th = stats.norm.ppf(0.90, loc=mu, scale=sigma)

# Z-score for a score of 90
z = (90 - mu) / sigma

print(f"P(X <= 70):        {prob_below_70:.4f}")
print(f"90th percentile:   {score_90th:.1f}")
print(f"Z-score for 90:    {z:.2f}")

P(X <= 70):        0.3556
90th percentile:   91.8
Z-score for 90:    1.14

Binomial Distribution

The binomial distribution models the number of successes in n independent trials, each with probability p. It requires fixed n, binary outcomes, constant p, and independence.

# A tutored student has a 70% chance of passing each quiz.
# In 8 quizzes, what is the probability of passing exactly 6?
n, p = 8, 0.70

prob_exactly_6 = stats.binom.pmf(6, n, p)
prob_at_most_5 = stats.binom.cdf(5, n, p)
prob_at_least_6 = 1 - stats.binom.cdf(5, n, p)

print(f"P(X = 6):   {prob_exactly_6:.4f}")
print(f"P(X <= 5):  {prob_at_most_5:.4f}")
print(f"P(X >= 6):  {prob_at_least_6:.4f}")

P(X = 6):   0.2965
P(X <= 5):  0.4482
P(X >= 6):  0.5518

Sampling and Statistical Inference

Standard Error

The standard error of the mean (SE) measures how much sample means vary from sample to sample. It decreases as sample size grows, because larger samples give more stable estimates.

n  = len(scores)
s  = scores.std()      # sample SD, ddof=1
se = s / np.sqrt(n)

print(f"n = {n}, s = {s:.2f}, SE = {se:.2f}")
# n = 10, s = 13.21, SE = 4.18

Confidence Intervals

A confidence interval gives a range of plausible values for a population parameter based on your sample. A 95% confidence interval means that if you repeated the study many times, 95% of the intervals constructed that way would contain the true population parameter. It does not mean there is a 95% probability that this specific interval contains the true value — once the interval is computed, the parameter is either in it or not.

# scipy.stats.t.interval gives the CI directly
confidence_level = 0.95
df = n - 1

ci_low, ci_high = stats.t.interval(
    confidence=confidence_level,
    df=df,
    loc=scores.mean(),
    scale=stats.sem(scores)      # scipy sem uses ddof=1
)

print(f"95% CI: ({ci_low:.2f}, {ci_high:.2f})")
# 95% CI: (65.45, 84.35)

Interpreted: based on this sample, a plausible range for the true population mean exam score is approximately 65.5 to 84.4. The width reflects the sample size (n=10) and variability (s=13.21). A larger sample would produce a narrower interval.

Hypothesis Testing in Python

Hypothesis testing uses sample data to decide whether there is enough evidence to reject a specific claim about a population. Every test has the same structure: a null hypothesis (H₀), an alternative hypothesis (H₁), a test statistic, a p-value, and a decision.

One-Sample t-Test

Use a one-sample t-test when you want to compare a sample mean to a known or hypothesized population value. Learn more in the dedicated one-sample t-test guide.

t_stat, p_value = stats.ttest_1samp(scores, popmean=70)
print(f"t = {t_stat:.3f}, p = {p_value:.4f}")
# t = 1.171, p = 0.2718

Independent-Samples t-Test (Welch)

When comparing means from two separate groups, use an independent-samples t-test. Welch's version (the default in SciPy) does not assume the groups have equal variance, which makes it the safer default. Full details in the two-sample t-test guide.

control = data[data["group"] == "Control"]["exam_score"]
tutored = data[data["group"] == "Tutored"]["exam_score"]

t_stat, p_value = stats.ttest_ind(control, tutored, equal_var=False)
print(f"Control mean: {control.mean()}, Tutored mean: {tutored.mean()}")
print(f"t = {t_stat:.3f}, p = {p_value:.4f}")
# Control mean: 66.2, Tutored mean: 83.6
# t = -3.041, p = 0.0163

Welch t(7.5) = -3.041, p = 0.016.
Since p < 0.05, we reject H0. The tutored group scored significantly
higher than the control group (83.6 vs 66.2 points, on average).

Paired t-Test

When each observation in group 1 is matched to a specific observation in group 2 — pre/post measurements on the same person, matched pairs — use the paired t-test.

# Suppose these are pre-test and post-test scores for the same 5 students
pre_test  = np.array([60, 65, 70, 55, 72])
post_test = np.array([68, 71, 75, 63, 79])

t_stat, p_value = stats.ttest_rel(pre_test, post_test)
mean_diff = (post_test - pre_test).mean()

print(f"Mean improvement: {mean_diff:.1f} points")
print(f"t = {t_stat:.3f}, p = {p_value:.4f}")
# Mean improvement: 7.2 points
# t = -8.367, p = 0.0011

Chi-Square Test

The chi-square test of independence checks whether two categorical variables are associated.

# Does pass/fail rate differ between Control and Tutored groups?
# Pass = score >= 70
data["pass_fail"] = (data["exam_score"] >= 70).map({True: "Pass", False: "Fail"})

contingency = pd.crosstab(data["group"], data["pass_fail"])
print(contingency)

chi2, p_val, dof, expected = stats.chi2_contingency(contingency)
print(f"χ² = {chi2:.3f}, p = {p_val:.4f}, df = {dof}")

pass_fail  Fail  Pass
group
Control       3     2
Tutored       0     5

χ² = 3.333, p = 0.0679, df = 1

With n=10 the test has very limited power. Expected cell counts are below 5,
which violates chi-square assumptions — interpret cautiously.

One-Way ANOVA

When you have three or more independent groups and want to test whether any of the group means differ, use one-way ANOVA. The F-test tells you whether at least one group mean is different, but not which ones. For that, you need a post-hoc test.

# Adding a third group (Online) for this example
online = np.array([69, 73, 71, 75, 70])
control_arr = control.values

f_stat, p_value = stats.f_oneway(control_arr, tutored.values, online)
print(f"F = {f_stat:.3f}, p = {p_value:.4f}")
# F = 7.981, p = 0.0051

Nonparametric Tests

Nonparametric tests make fewer distributional assumptions. They are not simply "tests for non-normal data" — each test has its own target and assumptions, often framed around medians or rank-based comparisons rather than means.

Correlation in Python

Correlation measures the direction and strength of the linear relationship between two variables. It ranges from −1 (perfect negative linear relationship) to +1 (perfect positive linear relationship). Zero indicates no linear relationship, though there may still be a nonlinear one. Learn more about Pearson correlation in the main topic guide.

x = data["hours_studied"]
y = data["exam_score"]

# Pearson: measures linear association
r_pearson, p_pearson = stats.pearsonr(x, y)

# Spearman: measures monotonic association (rank-based)
r_spearman, p_spearman = stats.spearmanr(x, y)

# Kendall tau: another rank-based measure, more robust with small samples
tau, p_kendall = stats.kendalltau(x, y)

print(f"Pearson r  = {r_pearson:.3f}, p = {p_pearson:.4f}")
print(f"Spearman r = {r_spearman:.3f}, p = {p_spearman:.4f}")
print(f"Kendall τ = {tau:.3f}, p = {p_kendall:.4f}")

Pearson r  = 0.975, p = 0.0000
Spearman r = 0.976, p = 0.0000
Kendall τ = 0.911, p = 0.0001

All three measures agree: there is a very strong positive association between hours studied and exam scores in this dataset. Use Pearson when both variables are continuous and the relationship appears linear. Use Spearman or Kendall when data are ordinal, the relationship is monotonic but not necessarily linear, or outliers are a concern.

Regression in Python

Simple Linear Regression

Simple linear regression models the relationship between one predictor (X) and a continuous outcome (Y) as a straight line.

model = smf.ols("exam_score ~ hours_studied", data=data).fit()
print(model.summary())

                 coef    std err       t      P>|t|   [0.025   0.975]
Intercept      34.0882     4.321    7.889    0.000    24.12    44.06
hours_studied   6.4706     0.659    9.823    0.000     4.95     7.99

R-squared:       0.9233
Adj. R-squared:  0.9137
F-statistic:     96.49  (p = 4.56e-06)

Reading this output:

• Intercept (34.09): the predicted exam score when hours studied is zero. In context this is not very meaningful, but it anchors the line.
• hours_studied coefficient (6.47): for each additional hour studied, the predicted exam score increases by 6.47 points, holding nothing else constant (since there is only one predictor here).
• p-values (both 0.000): both coefficients are statistically significant at α = 0.05.
• R² = 0.923: 92.3% of the variation in exam scores is explained by hours studied in this model. Explore R-squared in depth.
• 95% CI for slope (4.95, 7.99): a plausible range for the true slope. Excludes zero, consistent with the significant p-value.

Make predictions and examine residuals:

predicted = model.fittedvalues
residuals = model.resid

# Predict for a student who studied 6.5 hours
new_data = pd.DataFrame({"hours_studied": [6.5]})
prediction = model.predict(new_data)
print(f"Predicted score for 6.5 hours: {prediction[0]:.1f}")
# Predicted score for 6.5 hours: 76.1

Multiple Linear Regression

When you have more than one predictor, you fit a multiple regression model. Each coefficient represents the expected change in Y for a one-unit increase in that predictor, holding all other predictors in the model constant.

# C(group) tells statsmodels to treat "group" as a categorical variable
# The "Control" group becomes the reference category
mlr = smf.ols("exam_score ~ hours_studied + C(group)", data=data).fit()
print(mlr.summary())

In the output, the coefficient for C(group)[T.Tutored] estimates the average score difference between Tutored and Control students after accounting for study hours. This is the partial effect of tutoring — the adjusted difference once you have statistically controlled for how much each student studied. See our guide to multiple linear regression for more on interpretation, adjusted R², VIF, and diagnostics.

Logistic Regression

When the outcome is binary (pass/fail, yes/no), use logistic regression. It models the log-odds of the outcome as a linear function of the predictors.

# Create binary outcome: 1 = passed (score >= 70), 0 = failed
data["passed"] = (data["exam_score"] >= 70).astype(int)

logit_model = smf.logit("passed ~ hours_studied", data=data).fit()
print(logit_model.summary())

# Odds ratios (exponentiate the coefficients)
odds_ratios = np.exp(logit_model.params)
print(f"\nOdds ratio for hours_studied: {odds_ratios['hours_studied']:.3f}")

An odds ratio greater than 1 for hours_studied means that each additional hour studied is associated with higher odds of passing. Do not interpret logistic regression coefficients as linear changes in probability — convert to odds ratios or use predict() to get predicted probabilities.

Residual Diagnostics in Python

Fitting a model is only half the work. Checking residuals tells you whether the model's assumptions are reasonable and whether any observations are pulling the results disproportionately. These concepts are covered in depth in the statistical assumptions section.

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Residuals vs fitted values — look for patterns
axes[0].scatter(model.fittedvalues, model.resid, alpha=0.7)
axes[0].axhline(0, color="red", linestyle="--")
axes[0].set_xlabel("Fitted values")
axes[0].set_ylabel("Residuals")
axes[0].set_title("Residuals vs Fitted")

# Q-Q plot — residuals vs normal distribution quantiles
sm.qqplot(model.resid, line="s", ax=axes[1], alpha=0.7)
axes[1].set_title("Q-Q Plot of Residuals")

plt.tight_layout()
plt.show()

In the residuals-vs-fitted plot, look for points scattered randomly around zero with no obvious pattern. A curve or funnel shape suggests a nonlinear relationship or heteroscedasticity. In the Q-Q plot, points near the diagonal line suggest the residuals are approximately normally distributed, which matters for the validity of inference in small samples. See the Q-Q plots guide for more detail.

Statistical Visualization in Python

Good visualizations surface patterns that summary statistics can miss. Use Seaborn for statistical plots with sensible defaults, and Matplotlib when you need precise control. The data visualization section covers when to use each chart type.

Histogram

A histogram shows the distribution of a continuous variable: center, spread, skewness, and whether the data are approximately unimodal. The KDE line is a smoothed density estimate.

sns.histplot(data["exam_score"], kde=True, bins=6, color="steelblue")
plt.xlabel("Exam Score")
plt.ylabel("Count")
plt.title("Distribution of Exam Scores")
plt.show()

Box Plot

A box plot shows Q1, median, Q3, IQR, whiskers (typically 1.5 IQR), and potential outliers. It is especially useful for comparing groups side by side.

sns.boxplot(data=data, x="group", y="exam_score", palette="Blues")
plt.xlabel("Study Group")
plt.ylabel("Exam Score")
plt.title("Exam Score Distribution by Group")
plt.show()

Scatter Plot

A scatter plot shows the relationship between two continuous variables. Look for direction (positive/negative), form (linear/curved), strength, clusters, and outliers.

sns.regplot(data=data, x="hours_studied", y="exam_score",
            scatter_kws={"alpha":0.7}, line_kws={"color":"red"})
plt.xlabel("Hours Studied")
plt.ylabel("Exam Score")
plt.title("Hours Studied vs Exam Score")
plt.show()

SciPy vs statsmodels vs scikit-learn

A question that comes up constantly: which library should I use for this statistical analysis? They are not interchangeable.

A standard academic or business analysis uses SciPy for hypothesis tests and statsmodels for regression. A machine learning pipeline uses scikit-learn for training and evaluation. Many data science projects use all three.

End-to-End Analysis: Does Study Time Predict Exam Performance?

This section walks through a complete statistical workflow on the dataset defined earlier, showing how each step connects to the next.

Interactive Descriptive Statistics Calculator

Enter a list of numbers separated by commas to calculate descriptive statistics. The calculator uses sample formulas (ddof=1) for standard deviation and variance, matching the pandas default and the convention for analyzing samples rather than full populations.

Common Python Statistics Mistakes

Statistical Test Selection Guide

Method selection depends on your research question, measurement scale, study design, and assumptions — not on a single rule like "if normal, use t-test."

Effect Sizes in Python

Statistical significance tells you whether an effect exists in the sample. Effect size tells you how large it is. These are separate questions, and both matter. A study with n=10,000 can find statistically significant but trivially small effects.

# Cohen's d = (mean1 - mean2) / pooled SD
def cohens_d(g1, g2):
    n1, n2 = len(g1), len(g2)
    pooled_sd = np.sqrt(
        ((n1 - 1) * np.var(g1, ddof=1) + (n2 - 1) * np.var(g2, ddof=1))
        / (n1 + n2 - 2)
    )
    return (np.mean(g1) - np.mean(g2)) / pooled_sd

d = cohens_d(tutored.values, control.values)
print(f"Cohen's d = {d:.3f}")
# Cohen's d = 2.720  (very large effect)

Cohen's d of approximately 0.2 is considered small, 0.5 medium, and 0.8 large by conventional benchmarks, though these benchmarks vary by field. A d of 2.72 indicates that the tutored group's mean is 2.72 pooled standard deviations above the control group's mean.

How to Learn Statistics With Python

A practical learning sequence — each step builds on the previous one.

Python vs R for Statistics

This comes down to workflow, not capability. Both environments can perform virtually every classical statistical analysis.

Many researchers and data scientists work fluently in both. If your primary work is machine learning or software engineering alongside statistics, Python is the natural center. If your primary work is statistical modeling and publication-quality analysis, either works well and R may have advantages for specialized methods.

Glossary

Frequently Asked Questions

Sources and Further Reading