Intermediate

NumPy handles arrays. Pandas handles tables. But when you need to solve a system of equations, find the minimum of a complex function, integrate a curve, or run a statistical test, you need SciPy. It is the Swiss Army knife of scientific computing in Python, built on top of NumPy and packed with algorithms that scientists, engineers, and data professionals use every day.

SciPy is organised into subpackages, each covering a different domain: optimization, linear algebra, statistics, signal processing, interpolation, and more. You rarely import all of SciPy at once — instead, you pull in just the subpackage you need, keeping your code clean and your imports explicit.

In this tutorial, you will learn how to install SciPy, solve optimization problems, perform statistical tests, work with linear algebra, integrate functions numerically, interpolate data points, and process signals. By the end, you will have a practical toolkit for tackling real scientific and engineering problems in Python.

Written by Pubs

Python developer and educator with 15+ years building production systems across data engineering, web APIs, and AI tooling. Founder of Python How To Program — 270+ in-depth tutorials covering the modern Python stack.

View all tutorials by Pubs →

SciPy for Scientific Computing: Quick Example

Let us start with a common task: finding the minimum of a mathematical function. This comes up everywhere from machine learning (gradient descent) to engineering (minimising material cost).

# quick_scipy.py
from scipy.optimize import minimize
import numpy as np

# Define a function: the Rosenbrock function (a classic test problem)
def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Find the minimum starting from an initial guess
result = minimize(rosenbrock, x0=[0, 0], method='Nelder-Mead')

print(f"Minimum found at: x={result.x[0]:.4f}, y={result.x[1]:.4f}")
print(f"Minimum value: {result.fun:.6f}")
print(f"Converged: {result.success}")
print(f"Iterations: {result.nit}")

Output:

Minimum found at: x=1.0000, y=1.0000
Minimum value: 0.000000
Converged: True
Iterations: 141

What this does step by step: The Rosenbrock function is a famous test problem with a global minimum at (1, 1). We pass it to scipy.optimize.minimize with an initial guess of (0, 0) and the Nelder-Mead algorithm (a gradient-free method). SciPy iterates 141 times and finds the exact minimum. This same function can minimise any callable Python function, making it invaluable for fitting models, calibrating parameters, and solving engineering problems.

Installing SciPy

SciPy depends on NumPy and ships with precompiled binaries for all major platforms, so installation is straightforward.

# Install SciPy
pip install scipy

# Verify installation and check version
python -c "import scipy; print(scipy.__version__)"

# SciPy subpackages are imported individually:
# from scipy import optimize, stats, linalg, integrate, interpolate, signal

Unlike some scientific libraries, SciPy installs cleanly on Windows, macOS, and Linux without needing a C compiler. The pip package includes optimised BLAS and LAPACK routines, so you get near-Fortran performance out of the box. If you are using Anaconda, SciPy comes preinstalled.

Optimization: Finding Minimums and Solving Equations

The scipy.optimize module is one of the most-used parts of SciPy. It handles curve fitting, root finding, and general optimization.

Curve Fitting with curve_fit

When you have experimental data and want to find the best parameters for a model, curve_fit is your go-to tool.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Generate noisy exponential decay data
np.random.seed(42)
x_data = np.linspace(0, 10, 50)
y_data = 3.5 * np.exp(-0.7 * x_data) + np.random.normal(0, 0.2, 50)

# Define the model function
def exponential_decay(x, amplitude, decay_rate):
    return amplitude * np.exp(-decay_rate * x)

# Fit the model to data
params, covariance = curve_fit(exponential_decay, x_data, y_data)
amplitude, decay_rate = params

print(f"Fitted amplitude: {amplitude:.3f} (true: 3.500)")
print(f"Fitted decay rate: {decay_rate:.3f} (true: 0.700)")

# Calculate uncertainties from covariance matrix
uncertainties = np.sqrt(np.diag(covariance))
print(f"Amplitude uncertainty: +/- {uncertainties[0]:.3f}")
print(f"Decay rate uncertainty: +/- {uncertainties[1]:.3f}")

curve_fit returns two things: the optimal parameters and a covariance matrix. The diagonal of the covariance matrix gives you the variance of each parameter, and taking the square root gives you the standard error. This tells you not just the best fit, but how confident you should be in each parameter.

Root Finding

Finding where a function equals zero is fundamental to solving equations. SciPy offers several root-finding algorithms.

from scipy.optimize import brentq, fsolve
import numpy as np

# Find where x^3 - 2x - 5 = 0 using Brent's method (bracketed)
def cubic(x):
    return x**3 - 2*x - 5

root = brentq(cubic, 1, 3)  # Root must be between 1 and 3
print(f"Root of x^3 - 2x - 5: {root:.6f}")
print(f"Verification: f({root:.4f}) = {cubic(root):.2e}")

# Solve a system of nonlinear equations with fsolve
def system(variables):
    x, y = variables
    eq1 = x**2 + y**2 - 4    # Circle of radius 2
    eq2 = x - y**2 + 1        # Parabola
    return [eq1, eq2]

solution = fsolve(system, x0=[1, 1])
print(f"\nSystem solution: x={solution[0]:.4f}, y={solution[1]:.4f}")
print(f"Check circle: {solution[0]**2 + solution[1]**2:.4f} (should be 4)")

brentq is the fastest bracketed root finder — you give it an interval where the function changes sign, and it guarantees finding the root. fsolve handles systems of nonlinear equations, taking an initial guess and using Newton-type iterations to converge on the solution.

Statistical Analysis with scipy.stats

The scipy.stats module contains over 100 probability distributions plus a comprehensive collection of statistical tests.

from scipy import stats
import numpy as np

# Generate two samples
np.random.seed(42)
group_a = np.random.normal(loc=75, scale=10, size=100)  # Mean 75, std 10
group_b = np.random.normal(loc=78, scale=10, size=100)  # Mean 78, std 10

# Descriptive statistics
desc_a = stats.describe(group_a)
print(f"Group A: mean={desc_a.mean:.2f}, variance={desc_a.variance:.2f}")
print(f"  skewness={desc_a.skewness:.3f}, kurtosis={desc_a.kurtosis:.3f}")

# Independent samples t-test
t_stat, p_value = stats.ttest_ind(group_a, group_b)
print(f"\nT-test: t={t_stat:.3f}, p-value={p_value:.4f}")
print(f"Significant at 0.05? {'Yes' if p_value < 0.05 else 'No'}")

# Mann-Whitney U test (non-parametric alternative)
u_stat, p_mann = stats.mannwhitneyu(group_a, group_b, alternative='two-sided')
print(f"\nMann-Whitney U: U={u_stat:.1f}, p-value={p_mann:.4f}")

# Shapiro-Wilk normality test
w_stat, p_normal = stats.shapiro(group_a)
print(f"\nShapiro-Wilk (Group A): W={w_stat:.4f}, p-value={p_normal:.4f}")
print(f"Normal distribution? {'Yes' if p_normal > 0.05 else 'No'}")

# Pearson correlation
correlation, p_corr = stats.pearsonr(group_a[:50], group_b[:50])
print(f"\nPearson correlation: r={correlation:.3f}, p-value={p_corr:.4f}")

The t-test tells you whether two groups have significantly different means. The Mann-Whitney U test does the same thing but does not assume normal distributions. The Shapiro-Wilk test checks whether your data is normally distributed (important for choosing the right test). Always check your assumptions before picking a statistical test.

Working with Probability Distributions

from scipy import stats
import numpy as np

# Normal distribution
normal = stats.norm(loc=100, scale=15)  # IQ distribution: mean=100, std=15
print(f"P(IQ > 130): {1 - normal.cdf(130):.4f}")
print(f"P(85 < IQ < 115): {normal.cdf(115) - normal.cdf(85):.4f}")
print(f"IQ at 95th percentile: {normal.ppf(0.95):.1f}")

# Generate random samples
samples = normal.rvs(size=1000, random_state=42)
print(f"\nSample mean: {np.mean(samples):.1f}, std: {np.std(samples):.1f}")

# Fit a distribution to data
params = stats.norm.fit(samples)
print(f"Fitted params: mean={params[0]:.1f}, std={params[1]:.1f}")

# Chi-squared goodness of fit
observed = [18, 22, 20, 25, 15]
expected = [20, 20, 20, 20, 20]
chi2, p_chi = stats.chisquare(observed, expected)
print(f"\nChi-squared test: chi2={chi2:.2f}, p-value={p_chi:.4f}")

Every distribution in SciPy has the same interface: .pdf() for probability density, .cdf() for cumulative probability, .ppf() for percentiles (inverse CDF), .rvs() for random samples, and .fit() for parameter estimation. Once you learn one distribution, you know them all.

Linear Algebra with scipy.linalg

While NumPy has basic linear algebra, scipy.linalg adds decompositions, matrix functions, and specialised solvers that go well beyond the basics.

from scipy import linalg
import numpy as np

# Solve a system of linear equations: Ax = b
A = np.array([[3, 1, -1],
              [1, 4, 2],
              [2, 1, 3]])
b = np.array([4, 17, 13])

x = linalg.solve(A, b)
print(f"Solution: {x}")
print(f"Verification Ax: {A @ x}")

# LU decomposition
P, L, U = linalg.lu(A)
print(f"\nLU decomposition:")
print(f"L (lower triangular):\n{np.round(L, 3)}")
print(f"U (upper triangular):\n{np.round(U, 3)}")

# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = linalg.eig(A)
print(f"\nEigenvalues: {eigenvalues.real.round(3)}")

# Matrix determinant and inverse
det = linalg.det(A)
inv = linalg.inv(A)
print(f"\nDeterminant: {det:.1f}")
print(f"A @ inv(A) (should be identity):\n{np.round(A @ inv, 1)}")

# Singular Value Decomposition
U_svd, s, Vt = linalg.svd(A)
print(f"\nSingular values: {s.round(3)}")

linalg.solve is faster and more numerically stable than computing the inverse and multiplying. LU decomposition breaks a matrix into lower and upper triangular parts, which is useful when you need to solve the same system with different right-hand sides. SVD is the workhorse behind dimensionality reduction, recommendation systems, and image compression.

Numerical Integration

When you cannot find an analytical solution to an integral, numerical integration gives you an accurate approximation.

from scipy import integrate
import numpy as np

# Integrate a simple function: integral of sin(x) from 0 to pi
result, error = integrate.quad(lambda x: np.sin(x), 0, np.pi)
print(f"Integral of sin(x) from 0 to pi: {result:.6f} (exact: 2.0)")
print(f"Estimated error: {error:.2e}")

# Integrate a more complex function: Gaussian integral
result, error = integrate.quad(
    lambda x: np.exp(-x**2),
    -np.inf, np.inf  # Infinite limits are supported!
)
print(f"\nGaussian integral: {result:.6f} (exact: sqrt(pi) = {np.sqrt(np.pi):.6f})")

# Double integral: integral of x*y over the unit circle
def integrand(y, x):  # Note: y first, then x
    return x * y

def y_lower(x):
    return -np.sqrt(1 - x**2)

def y_upper(x):
    return np.sqrt(1 - x**2)

result, error = integrate.dblquad(integrand, -1, 1, y_lower, y_upper)
print(f"\nDouble integral of xy over unit circle: {result:.6f} (exact: 0)")

# Solve an ODE: dy/dt = -2y, y(0) = 1 (exponential decay)
def decay(t, y):
    return -2 * y

t_span = (0, 5)
t_eval = np.linspace(0, 5, 100)
solution = integrate.solve_ivp(decay, t_span, [1.0], t_eval=t_eval)

print(f"\ny(0) = {solution.y[0][0]:.4f}")
print(f"y(5) = {solution.y[0][-1]:.6f} (exact: {np.exp(-10):.6f})")

quad handles single integrals and even supports infinite limits. dblquad handles double integrals. solve_ivp solves initial value problems for ordinary differential equations, which is essential for modelling anything that changes over time: population growth, chemical reactions, mechanical systems, or circuit dynamics.

Interpolation: Filling in the Gaps

When you have data at discrete points and need values between them, interpolation creates a smooth function that passes through your data.

from scipy.interpolate import interp1d, CubicSpline
import numpy as np

# Known data points (e.g., temperature readings every 3 hours)
hours = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
temps = np.array([15, 13, 14, 18, 24, 26, 23, 19, 16])

# Linear interpolation
linear_interp = interp1d(hours, temps, kind='linear')

# Cubic spline interpolation (smoother)
cubic_interp = CubicSpline(hours, temps)

# Evaluate at every half hour
fine_hours = np.linspace(0, 24, 49)
linear_temps = linear_interp(fine_hours)
cubic_temps = cubic_interp(fine_hours)

# Compare at hour 7.5
print(f"Temperature at 7:30 AM:")
print(f"  Linear interpolation: {linear_interp(7.5):.1f} C")
print(f"  Cubic spline: {cubic_interp(7.5):.1f} C")

# Cubic spline also gives derivatives
print(f"\nRate of change at noon: {cubic_interp(12, 1):.2f} C/hour")
print(f"Rate of change at 6 PM: {cubic_interp(18, 1):.2f} C/hour")

Linear interpolation draws straight lines between points -- simple but produces sharp corners. Cubic splines create smooth curves that pass through every point and have continuous first and second derivatives. The CubicSpline object can also compute derivatives at any point, which is useful for finding rates of change.

Signal Processing Basics

The scipy.signal module handles filtering, spectral analysis, and signal manipulation -- essential for audio processing, sensor data, and time series analysis.

from scipy import signal
import numpy as np

# Create a noisy signal: clean sine wave + noise
np.random.seed(42)
fs = 1000  # Sampling frequency (Hz)
t = np.linspace(0, 1, fs, endpoint=False)
clean_signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 50 * t)
noisy_signal = clean_signal + 0.5 * np.random.randn(len(t))

# Design a low-pass Butterworth filter (remove frequencies above 20 Hz)
nyquist = fs / 2
cutoff = 20 / nyquist  # Normalize to Nyquist frequency
b, a = signal.butter(N=4, Wn=cutoff, btype='low')

# Apply the filter
filtered = signal.filtfilt(b, a, noisy_signal)

print(f"Original signal RMS: {np.sqrt(np.mean(noisy_signal**2)):.3f}")
print(f"Filtered signal RMS: {np.sqrt(np.mean(filtered**2)):.3f}")
print(f"Noise reduced by: {(1 - np.sqrt(np.mean((filtered - clean_signal[:len(filtered)])**2)) / np.sqrt(np.mean((noisy_signal - clean_signal)**2))) * 100:.1f}%")

# Find peaks in the filtered signal
peaks, properties = signal.find_peaks(filtered, height=0.5, distance=50)
print(f"\nPeaks found: {len(peaks)}")
print(f"Peak heights: {properties['peak_heights'][:5].round(3)}")

signal.butter designs a Butterworth filter with a smooth frequency response. signal.filtfilt applies it forward and backward to eliminate phase distortion. find_peaks locates local maxima in a signal, with parameters to control minimum height and distance between peaks. This pipeline -- design filter, apply filter, find features -- is the backbone of most signal processing workflows.

Real-World Example: Analysing Experimental Data

Let us combine multiple SciPy tools to analyse a realistic dataset: fitting a model to noisy measurements, testing hypotheses, and quantifying uncertainty.

import numpy as np
from scipy import stats, optimize, integrate

# Simulate experimental data: drug dosage vs response
np.random.seed(42)
doses = np.array([0, 0.5, 1, 2, 5, 10, 20, 50, 100])
true_response = 100 * doses / (10 + doses)  # Hill equation (pharmacology)
measured_response = true_response + np.random.normal(0, 5, len(doses))
measured_response = np.clip(measured_response, 0, 100)

# Fit the Hill equation to data
def hill_equation(dose, v_max, k_half):
    return v_max * dose / (k_half + dose)

params, cov = optimize.curve_fit(
    hill_equation, doses, measured_response,
    p0=[100, 10],  # Initial guesses
    bounds=([0, 0], [200, 100])  # Parameter bounds
)
v_max, k_half = params
errors = np.sqrt(np.diag(cov))

print("Hill Equation Fit Results:")
print(f"  V_max = {v_max:.1f} +/- {errors[0]:.1f} (true: 100)")
print(f"  K_half = {k_half:.1f} +/- {errors[1]:.1f} (true: 10)")

# Calculate R-squared
predicted = hill_equation(doses, *params)
ss_res = np.sum((measured_response - predicted)**2)
ss_tot = np.sum((measured_response - np.mean(measured_response))**2)
r_squared = 1 - ss_res / ss_tot
print(f"  R-squared = {r_squared:.4f}")

# Calculate Area Under the Curve (AUC) using integration
auc, auc_error = integrate.quad(
    lambda d: hill_equation(d, *params), 0, 100
)
print(f"\nArea Under Curve (0-100): {auc:.1f} +/- {auc_error:.2e}")

# Statistical test: is the response at dose=50 significantly different from dose=5?
np.random.seed(42)
samples_dose5 = hill_equation(5, *params) + np.random.normal(0, 5, 30)
samples_dose50 = hill_equation(50, *params) + np.random.normal(0, 5, 30)

t_stat, p_value = stats.ttest_ind(samples_dose5, samples_dose50)
print(f"\nDose 5 vs Dose 50 comparison:")
print(f"  Mean at dose 5: {np.mean(samples_dose5):.1f}")
print(f"  Mean at dose 50: {np.mean(samples_dose50):.1f}")
print(f"  t-statistic: {t_stat:.3f}")
print(f"  p-value: {p_value:.2e}")
print(f"  Significant difference: {'Yes' if p_value < 0.05 else 'No'}")

This example mirrors a real pharmacology workflow: fit a dose-response curve, quantify parameter uncertainty, calculate the area under the curve (a standard efficacy metric), and test whether two dosage levels produce significantly different responses. The same pattern applies to any field where you need to model data and draw statistical conclusions.

Frequently Asked Questions

What is the difference between NumPy and SciPy?

NumPy provides the fundamental array data structure and basic operations like element-wise math, reshaping, and basic linear algebra. SciPy builds on NumPy and adds higher-level scientific algorithms: optimization, statistics, signal processing, integration, interpolation, and advanced linear algebra. Think of NumPy as the foundation and SciPy as the specialised toolkit.

Should I use scipy.linalg or numpy.linalg?

scipy.linalg is a superset of numpy.linalg with additional decompositions and solvers. It also always uses BLAS and LAPACK, which can be faster. For basic operations like dot or norm, either works fine. For decompositions (LU, Cholesky, SVD) or specialised solvers, prefer scipy.linalg.

How do I choose between different optimization methods?

If your function is smooth and differentiable, use L-BFGS-B or BFGS (fast, gradient-based). If you cannot compute gradients, use Nelder-Mead or Powell. If your function has many local minima, consider differential_evolution (global optimizer). For bounded problems, L-BFGS-B handles box constraints natively.

Can SciPy handle large datasets?

SciPy works with NumPy arrays, so it handles arrays that fit in memory efficiently. For very large sparse matrices, use scipy.sparse, which stores only non-zero elements. For datasets larger than memory, consider chunked processing or libraries like Dask that parallelize SciPy operations.

How do I choose the right statistical test?

For comparing two group means with normal data, use the t-test (ttest_ind). For non-normal data, use Mann-Whitney U (mannwhitneyu). For more than two groups, use one-way ANOVA (f_oneway) or Kruskal-Wallis (kruskal). Always check normality with shapiro first and check equal variances with levene before choosing a parametric test.

Wrapping Up

SciPy gives you access to decades of scientific computing algorithms through a clean, consistent Python interface. You have learned how to optimise functions and fit models with scipy.optimize, run statistical tests and work with distributions using scipy.stats, solve linear algebra problems with scipy.linalg, integrate functions numerically with scipy.integrate, interpolate data with scipy.interpolate, and process signals with scipy.signal. The key to using SciPy effectively is knowing which subpackage to reach for and understanding the assumptions behind each algorithm. Start with the examples in this tutorial, adapt them to your own data, and you will find that SciPy handles the mathematical heavy lifting while you focus on the science.

Related Articles

• How To Use NumPy Arrays in Python
• How To Use Pandas for Data Analysis in Python
• How To Create Data Visualizations with Seaborn in Python

Intermediate

Raw numbers in a spreadsheet rarely tell a compelling story. A well-crafted chart, on the other hand, can reveal patterns in seconds that would take minutes of scanning rows and columns. Python’s seaborn library sits on top of matplotlib and turns complex statistical visualizations into one-line function calls with beautiful default styling. Whether you need a quick histogram, a correlation heatmap, or a multi-faceted regression plot, Seaborn handles the heavy lifting so you can focus on understanding your data.

In this tutorial, you will learn how to install Seaborn, create common chart types including scatter plots, bar charts, histograms, and heatmaps, customise styles and colour palettes, work with real-world datasets, build multi-plot grids with FacetGrid, and export publication-ready figures. By the end, you will have the tools to turn any pandas DataFrame into a visual story.

Written by Pubs

Python developer and educator with 15+ years building production systems across data engineering, web APIs, and AI tooling. Founder of Python How To Program — 270+ in-depth tutorials covering the modern Python stack.

View all tutorials by Pubs →

Data Visualizations with Seaborn: Quick Example

Let us start with a scatter plot that reveals the relationship between two variables in a built-in dataset. This gets you from zero to a polished chart in four lines.

# quick_seaborn.py
import seaborn as sns
import matplotlib.pyplot as plt

tips = sns.load_dataset("tips")
sns.scatterplot(data=tips, x="total_bill", y="tip", hue="time")
plt.title("Tips by Total Bill")
plt.tight_layout()
plt.savefig("tips_scatter.png", dpi=150)
plt.show()

What this does step by step:

sns.load_dataset("tips") fetches a built-in DataFrame with restaurant tipping data. sns.scatterplot creates a scatter plot with total_bill on the x-axis and tip on the y-axis, automatically colouring points by time (Lunch vs Dinner). The plt.tight_layout() call prevents labels from being clipped, and savefig exports the chart at 150 DPI. You should see a clear upward trend: larger bills tend to produce larger tips, with dinner and lunch points neatly separated by colour.

Installing Seaborn and Its Dependencies

Seaborn requires matplotlib, pandas, and numpy, but pip handles all of them automatically.

# Install seaborn (pulls matplotlib, pandas, numpy)
pip install seaborn

# Verify the installation
python -c "import seaborn as sns; print(sns.__version__)"

If you are using Jupyter notebooks, Seaborn works out of the box. For scripts, remember to call plt.show() to display charts in a window, or use plt.savefig() to save them directly to files. Seaborn version 0.12 and above uses a new interface with the objects module, but the classic functional API we use here remains fully supported and is still the most common approach you will find in tutorials and production code.

Essential Chart Types in Seaborn

Seaborn organises its plotting functions into categories based on the type of relationship you want to show. Understanding which function to use for your data is half the battle.

Scatter Plots and Line Plots for Relationships

Use sns.scatterplot when you want to see how two continuous variables relate. Use sns.lineplot when your x-axis represents a sequence like time.

import seaborn as sns
import matplotlib.pyplot as plt

# Scatter plot with size encoding
tips = sns.load_dataset("tips")
sns.scatterplot(
    data=tips,
    x="total_bill",
    y="tip",
    hue="day",
    size="size",
    sizes=(20, 200),
    alpha=0.7
)
plt.title("Tip Amount vs Total Bill by Day")
plt.show()

# Line plot for time series
fmri = sns.load_dataset("fmri")
sns.lineplot(
    data=fmri,
    x="timepoint",
    y="signal",
    hue="region",
    style="event",
    errorbar="sd"
)
plt.title("FMRI Signal Over Time")
plt.show()

The hue parameter assigns colours by category, size maps a numeric column to marker diameter, and alpha controls transparency so overlapping points remain visible. For line plots, errorbar="sd" adds a shaded band showing the standard deviation, giving you a sense of how much the data varies at each time point.

Bar Charts and Count Plots for Categories

When one of your axes represents a category rather than a continuous number, bar charts are the right choice.

# Average tip by day
sns.barplot(data=tips, x="day", y="tip", hue="sex", errorbar="sd")
plt.title("Average Tip by Day and Gender")
plt.show()

# Count occurrences in each category
sns.countplot(data=tips, x="day", hue="smoker", palette="Set2")
plt.title("Visits per Day by Smoker Status")
plt.show()

sns.barplot automatically computes the mean and adds error bars. sns.countplot simply tallies how many rows fall into each category, which is perfect for understanding the distribution of categorical variables in your dataset.

Histograms and Distribution Plots

Understanding how a single variable is distributed is fundamental to any analysis. Seaborn gives you several options.

# Histogram with KDE overlay
sns.histplot(data=tips, x="total_bill", bins=25, kde=True, color="steelblue")
plt.title("Distribution of Total Bills")
plt.show()

# KDE plot comparing groups
sns.kdeplot(data=tips, x="total_bill", hue="time", fill=True, alpha=0.5)
plt.title("Bill Distribution: Lunch vs Dinner")
plt.show()

# Box plot for comparing distributions
sns.boxplot(data=tips, x="day", y="total_bill", hue="smoker", palette="coolwarm")
plt.title("Bill Amounts by Day and Smoker Status")
plt.show()

Setting kde=True overlays a kernel density estimate curve on the histogram, smoothing out the bars into a continuous shape. sns.kdeplot with fill=True creates shaded density curves, making it easy to compare two groups visually. Box plots show the median, quartiles, and outliers in a compact format that works well when you have multiple categories to compare side by side.

Building Heatmaps for Correlation Analysis

Heatmaps turn a matrix of numbers into a colour-coded grid, making it easy to spot strong correlations at a glance. They are one of Seaborn’s most popular features.

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# Load a dataset and compute correlation matrix
penguins = sns.load_dataset("penguins").dropna()
numeric_cols = penguins.select_dtypes(include=[np.number])
corr_matrix = numeric_cols.corr()

# Create an annotated heatmap
plt.figure(figsize=(8, 6))
sns.heatmap(
    corr_matrix,
    annot=True,
    fmt=".2f",
    cmap="RdBu_r",
    center=0,
    square=True,
    linewidths=0.5,
    vmin=-1,
    vmax=1
)
plt.title("Penguin Measurements Correlation Matrix")
plt.tight_layout()
plt.show()

The annot=True parameter prints the correlation value inside each cell, while fmt=".2f" rounds it to two decimal places. cmap="RdBu_r" uses a red-blue diverging colour scheme where red means strong positive correlation and blue means strong negative. Setting center=0 ensures that zero correlation appears as white, making the pattern immediately interpretable.

Customising Styles and Colour Palettes

Seaborn comes with five built-in themes and a flexible palette system that lets you match your charts to any brand or presentation style.

# Set a global theme
sns.set_theme(style="whitegrid", font_scale=1.2)

# Compare available styles
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
styles = ["darkgrid", "whitegrid", "dark", "white"]

for ax, style in zip(axes, styles):
    with sns.axes_style(style):
        sns.barplot(data=tips, x="day", y="tip", ax=ax, palette="viridis")
        ax.set_title(style)

plt.tight_layout()
plt.show()

The five styles are darkgrid, whitegrid, dark, white, and ticks. For colour palettes, you can use named palettes like "viridis", "Set2", or "coolwarm", or create your own with sns.color_palette("husl", 8) for 8 evenly spaced hues. The font_scale parameter is especially useful when preparing charts for presentations where you need larger text.

# Custom colour palette
custom_palette = sns.color_palette(["#2ecc71", "#e74c3c", "#3498db", "#f39c12"])
sns.barplot(data=tips, x="day", y="tip", palette=custom_palette)
plt.title("Tips by Day (Custom Colours)")
plt.show()

Multi-Plot Grids with FacetGrid

When you want to see how a pattern changes across different categories, FacetGrid creates a grid of small multiples, each showing the same chart type for a different subset of your data.

# FacetGrid: histogram for each day
g = sns.FacetGrid(tips, col="day", col_wrap=2, height=4)
g.map_dataframe(sns.histplot, x="total_bill", kde=True, color="steelblue")
g.set_titles("{col_name}")
g.set_axis_labels("Total Bill ($)", "Count")
g.tight_layout()
plt.show()

# PairGrid: scatter matrix for numeric columns
penguins = sns.load_dataset("penguins").dropna()
g = sns.pairplot(
    penguins,
    hue="species",
    diag_kind="kde",
    plot_kws={"alpha": 0.6, "s": 40}
)
g.fig.suptitle("Penguin Species Comparison", y=1.02)
plt.show()

FacetGrid takes a DataFrame and a column to split on (col for columns, row for rows). The col_wrap parameter controls how many plots fit in each row before wrapping. pairplot is a convenience function that creates a scatter matrix showing every numeric column against every other, with distribution plots on the diagonal. It is one of the fastest ways to explore a new dataset.

Real-World Example: Analysing Flight Delays

Let us put everything together with a real-world scenario. Suppose you have a dataset of flight information and want to understand seasonal passenger patterns.

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

# Load the flights dataset (monthly passenger counts 1949-1960)
flights = sns.load_dataset("flights")

# Pivot for heatmap
flights_pivot = flights.pivot(index="month", columns="year", values="passengers")

# Create a comprehensive dashboard
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Heatmap of passengers by month and year
sns.heatmap(flights_pivot, annot=True, fmt="d", cmap="YlOrRd",
            ax=axes[0, 0], cbar_kws={"label": "Passengers"})
axes[0, 0].set_title("Monthly Passengers (1949-1960)")

# 2. Line plot showing yearly trends
sns.lineplot(data=flights, x="year", y="passengers", hue="month",
             palette="tab20", ax=axes[0, 1], legend=False)
axes[0, 1].set_title("Passenger Trends by Month")

# 3. Box plot of monthly distributions
sns.boxplot(data=flights, x="month", y="passengers",
            palette="coolwarm", ax=axes[1, 0])
axes[1, 0].set_title("Monthly Passenger Distribution")
axes[1, 0].tick_params(axis="x", rotation=45)

# 4. Bar plot of yearly totals
yearly = flights.groupby("year")["passengers"].sum().reset_index()
sns.barplot(data=yearly, x="year", y="passengers",
            palette="Blues_d", ax=axes[1, 1])
axes[1, 1].set_title("Total Passengers per Year")
axes[1, 1].tick_params(axis="x", rotation=45)

plt.suptitle("Flight Passenger Analysis Dashboard", fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig("flight_dashboard.png", dpi=150, bbox_inches="tight")
plt.show()

This dashboard combines four chart types into a single figure. The heatmap reveals that summer months consistently have the highest passenger counts. The line plot shows a clear upward trend across all months over the years. The box plot highlights that July has the widest range of values, while the bar chart confirms that total yearly passengers grew steadily from 1949 to 1960. Building multi-panel dashboards like this is one of the most practical skills you can develop with Seaborn.

Saving and Exporting Publication-Ready Figures

Creating a great chart is only half the job. You also need to export it at the right resolution and format for your audience.

# Save as PNG at high resolution
fig, ax = plt.subplots(figsize=(10, 6))
sns.histplot(data=tips, x="total_bill", kde=True, ax=ax)
ax.set_title("Distribution of Total Bills")

# PNG for web and presentations
plt.savefig("chart.png", dpi=300, bbox_inches="tight", facecolor="white")

# SVG for scalable vector graphics (papers, reports)
plt.savefig("chart.svg", format="svg", bbox_inches="tight")

# PDF for LaTeX documents
plt.savefig("chart.pdf", format="pdf", bbox_inches="tight")

plt.close()
print("Charts saved successfully")

Use dpi=300 for print quality and dpi=150 for web use. The bbox_inches="tight" parameter trims whitespace around the chart. SVG format is ideal for reports because it scales without pixelation. Always call plt.close() after saving to free memory, especially when generating many charts in a loop.

Frequently Asked Questions

What is the difference between Seaborn and Matplotlib?

Matplotlib is the foundation that handles the actual drawing of axes, lines, and shapes. Seaborn sits on top of matplotlib and provides a higher-level interface with better default styles, built-in statistical aggregation, and simpler syntax for common chart types. You can mix both in the same script since every Seaborn function returns a matplotlib axes object.

Can I use Seaborn with data that is not in a pandas DataFrame?

Yes, most Seaborn functions accept numpy arrays, Python lists, or dictionaries in addition to DataFrames. However, the DataFrame interface is the most powerful because it allows you to reference column names directly for parameters like hue, size, and style. Converting your data to a DataFrame first is almost always worth the extra line of code.

How do I change the figure size in Seaborn?

For axes-level functions like sns.scatterplot, create the figure first with plt.figure(figsize=(10, 6)) or fig, ax = plt.subplots(figsize=(10, 6)) and pass the ax parameter. For figure-level functions like sns.catplot, use the height and aspect parameters directly.

Why do my Seaborn charts look different from the examples online?

Seaborn version 0.12 changed several default behaviours, including the default theme and some function names. Run sns.set_theme() at the start of your script to apply the modern defaults, and check your version with sns.__version__. Also note that some older tutorials use deprecated functions like distplot which was replaced by histplot and kdeplot.

How do I add labels and annotations to Seaborn plots?

Since Seaborn returns matplotlib axes, you can use all matplotlib annotation functions. Call ax.set_xlabel(), ax.set_ylabel(), and ax.set_title() for basic labels. For annotations pointing to specific data points, use ax.annotate("text", xy=(x, y), xytext=(x2, y2), arrowprops=dict(arrowstyle="->")).

Wrapping Up

Seaborn transforms the often tedious process of data visualization into something fast and enjoyable. You have learned how to create scatter plots, bar charts, histograms, heatmaps, and multi-plot grids, all with Seaborn’s clean one-line syntax. The key to mastering Seaborn is understanding which chart type matches your data: use scatter plots for two continuous variables, bar charts for categorical comparisons, histograms and KDE plots for distributions, and heatmaps for correlation matrices. Combined with FacetGrid for multi-panel layouts and Seaborn’s built-in themes for consistent styling, you now have a complete toolkit for turning raw data into compelling visual stories. Start with the built-in datasets to practice, then apply these patterns to your own data.

Related Articles

• How To Use Pandas for Data Analysis in Python
• How To Use NumPy Arrays in Python
• How To Read and Write CSV Files in Python

Intermediate

You have probably written hundreds of Python functions, but have you ever wondered what happens when a function is defined inside another function — and the inner function remembers the outer function’s variables even after the outer function has finished running? That is a closure, and it is one of the most powerful and underused features in Python.

Closures let you create lightweight, stateful functions without defining a class. They are used extensively in decorators, callback systems, event handlers, and factory patterns. Understanding closures also unlocks a deeper understanding of how Python’s scoping rules actually work.

In this tutorial, you will learn how nested functions work, what closures are and how Python creates them, the LEGB scoping rule, the nonlocal keyword for modifying enclosed variables, and practical patterns where closures replace classes. By the end, you will be using closures confidently in your own projects.

Written by Pubs

Python developer and educator with 15+ years building production systems across data engineering, web APIs, and AI tooling. Founder of Python How To Program — 270+ in-depth tutorials covering the modern Python stack.

View all tutorials by Pubs →

Python Closures: Quick Example

Here is the simplest possible closure — a function that remembers a greeting prefix and uses it every time you call it.

# quick_closure.py
def make_greeter(prefix):
    def greet(name):
        return f"{prefix}, {name}!"
    return greet

hello = make_greeter("Hello")
howdy = make_greeter("Howdy")

print(hello("Alice"))
print(howdy("Bob"))
print(hello("Charlie"))

Output:

Hello, Alice!
Howdy, Bob!
Hello, Charlie!

The make_greeter function returns the inner greet function. Even though make_greeter has finished executing, the returned greet function still remembers the prefix value it was created with. That is a closure — the inner function “closes over” the variable from its enclosing scope. Let us explore how this works under the hood.

What Are Closures and Nested Functions?

A nested function (also called an inner function) is simply a function defined inside another function. The outer function is sometimes called the enclosing function. In Python, functions are first-class objects — you can pass them as arguments, return them from other functions, and assign them to variables.

A closure is a special case of a nested function: it is a function object that remembers values from its enclosing lexical scope even when the enclosing function is no longer active. For a closure to exist, three conditions must be met: there must be a nested function, the nested function must reference a variable from the enclosing scope, and the enclosing function must return the nested function.

You can verify that a function is a closure by inspecting its __closure__ attribute. If it returns None, the function is not a closure. If it returns a tuple of cell objects, each cell contains one of the free variables the function has closed over.

Understanding LEGB Scoping Rules

Python resolves variable names using the LEGB rule, checking each scope in order until it finds the variable. Understanding this rule is essential for understanding how closures capture variables.

# legb_demo.py
x = "global"

def outer():
    x = "enclosing"
    
    def inner():
        x = "local"
        print(f"Inner sees: {x}")
    
    inner()
    print(f"Outer sees: {x}")

outer()
print(f"Module sees: {x}")

Output:

Inner sees: local
Outer sees: enclosing
Module sees: global

LEGB stands for Local, Enclosing, Global, Built-in. When Python encounters a variable name, it first checks the local scope (inside the current function), then the enclosing scope (any outer functions), then the global scope (module level), and finally the built-in scope (Python’s built-in names like print and len). A closure captures variables from the Enclosing scope — the “E” in LEGB.

The nonlocal Keyword

By default, you can read enclosed variables from within a closure, but you cannot reassign them. If you try to assign a new value to an enclosed variable, Python creates a new local variable instead. The nonlocal keyword tells Python to look in the enclosing scope for the variable, allowing you to modify it.

# nonlocal_demo.py
def make_counter():
    count = 0
    
    def increment():
        nonlocal count
        count += 1
        return count
    
    def get_count():
        return count
    
    return increment, get_count

increment, get_count = make_counter()
print(increment())
print(increment())
print(increment())
print(f"Final count: {get_count()}")

Output:

1
2
3
Final count: 3

Without nonlocal count, the line count += 1 would raise an UnboundLocalError because Python would treat count as a local variable being referenced before assignment. The nonlocal declaration explicitly tells Python that count lives in the enclosing scope and should be modified there. This is what makes closures stateful — they can maintain and update state across calls.

Practical Closure Patterns

Closures are not just a theoretical concept — they solve real problems more elegantly than alternatives. Here are three patterns you will use regularly.

Factory Functions

A factory function creates and returns specialized functions. This is cleaner than creating a class when you just need a callable with some configuration baked in.

# factory_demo.py
def make_multiplier(factor):
    def multiply(number):
        return number * factor
    return multiply

double = make_multiplier(2)
triple = make_multiplier(3)
to_cents = make_multiplier(100)

print(f"Double 5: {double(5)}")
print(f"Triple 5: {triple(5)}")
print(f"$4.99 in cents: {to_cents(4.99)}")

# Apply to a list
prices = [9.99, 14.50, 3.25]
prices_in_cents = list(map(to_cents, prices))
print(f"Prices in cents: {prices_in_cents}")

Output:

Double 5: 10
Triple 5: 15
$4.99 in cents: 499.0
Prices in cents: [999.0, 1450.0, 325.0]

Memoization Cache

Closures can maintain a cache dictionary that persists across calls, implementing memoization without global variables or classes.

# memoize_demo.py
def memoize(func):
    cache = {}
    
    def wrapper(*args):
        if args not in cache:
            cache[args] = func(*args)
            print(f"  Computing {func.__name__}{args}")
        else:
            print(f"  Cache hit for {func.__name__}{args}")
        return cache[args]
    
    wrapper.cache = cache  # Expose cache for inspection
    return wrapper

@memoize
def fibonacci(n):
    if n < 2:
        return n
    return fibonacci(n - 1) + fibonacci(n - 2)

print(f"fib(6) = {fibonacci(6)}")
print(f"fib(4) = {fibonacci(4)}")
print(f"Cache size: {len(fibonacci.cache)} entries")

Output:

  Computing fibonacci(6,)
  Computing fibonacci(5,)
  Computing fibonacci(4,)
  Computing fibonacci(3,)
  Computing fibonacci(2,)
  Computing fibonacci(1,)
  Computing fibonacci(0,)
  Cache hit for fibonacci(1,)
  Cache hit for fibonacci(2,)
  Cache hit for fibonacci(3,)
  Cache hit for fibonacci(4,)
fib(6) = 8
  Cache hit for fibonacci(4,)
fib(4) = 3
Cache size: 7 entries

Event Handlers and Callbacks

Closures are perfect for callbacks that need context. Instead of creating a class just to hold one piece of state, create a closure.

# callback_demo.py
def make_logger(log_level):
    messages = []
    
    def log(message):
        entry = f"[{log_level.upper()}] {message}"
        messages.append(entry)
        print(entry)
    
    def get_logs():
        return messages.copy()
    
    def clear():
        nonlocal messages
        messages = []
    
    log.get_logs = get_logs
    log.clear = clear
    return log

error_log = make_logger("error")
info_log = make_logger("info")

error_log("Database connection failed")
error_log("Retry attempt 1")
info_log("Server started on port 8080")
error_log("Retry succeeded")

print(f"\nError log has {len(error_log.get_logs())} entries")
print(f"Info log has {len(info_log.get_logs())} entries")

Output:

[ERROR] Database connection failed
[ERROR] Retry attempt 1
[INFO] Server started on port 8080
[ERROR] Retry succeeded

Error log has 3 entries
Info log has 1 entries

Each logger maintains its own independent list of messages because each call to make_logger creates a new messages list in a new enclosing scope. This is the same isolation you would get from separate class instances, but with less boilerplate.

Closures vs Classes: When to Use Each

A common question is whether to use a closure or a class. Both can maintain state, but they have different strengths.

# closure_vs_class.py
# Closure approach
def make_accumulator_closure(initial=0):
    total = initial
    def add(value):
        nonlocal total
        total += value
        return total
    return add

# Class approach
class Accumulator:
    def __init__(self, initial=0):
        self.total = initial
    
    def add(self, value):
        self.total += value
        return self.total

# Both work the same way
closure_acc = make_accumulator_closure(10)
class_acc = Accumulator(10)

print(f"Closure: {closure_acc(5)}, {closure_acc(3)}")
print(f"Class:   {class_acc.add(5)}, {class_acc.add(3)}")

Output:

Closure: 15, 18
Class:   15, 18

The rule of thumb: if your "object" has one main action and minimal state, a closure is simpler. If it has multiple methods, complex state, or needs to participate in inheritance, use a class.

Real-Life Example: Building a Rate Limiter with Closures

Let us build a practical rate limiter that tracks function calls and enforces a maximum number of calls per time window. This demonstrates closures maintaining complex state across multiple calls.

# rate_limiter.py
import time

def rate_limit(max_calls, window_seconds):
    call_timestamps = []
    
    def decorator(func):
        def wrapper(*args, **kwargs):
            nonlocal call_timestamps
            now = time.time()
            # Remove timestamps outside the window
            call_timestamps = [t for t in call_timestamps if now - t < window_seconds]
            
            if len(call_timestamps) >= max_calls:
                wait_time = window_seconds - (now - call_timestamps[0])
                print(f"Rate limited! Try again in {wait_time:.1f}s")
                return None
            
            call_timestamps.append(now)
            remaining = max_calls - len(call_timestamps)
            print(f"Call allowed ({remaining} remaining in window)")
            return func(*args, **kwargs)
        
        wrapper.get_usage = lambda: len([t for t in call_timestamps if time.time() - t < window_seconds])
        return wrapper
    return decorator

@rate_limit(max_calls=3, window_seconds=5)
def fetch_data(query):
    return f"Results for: {query}"

# Simulate rapid API calls
for i in range(5):
    result = fetch_data(f"query_{i}")
    if result:
        print(f"  Got: {result}")
    time.sleep(0.5)

print(f"\nCurrent usage: {fetch_data.get_usage()} calls in window")

Output:

Call allowed (2 remaining in window)
  Got: Results for: query_0
Call allowed (1 remaining in window)
  Got: Results for: query_1
Call allowed (0 remaining in window)
  Got: Results for: query_2
Rate limited! Try again in 3.5s
Rate limited! Try again in 3.0s

Current usage: 3 calls in window

This rate limiter uses three levels of closures: rate_limit captures the configuration (max_calls, window_seconds), decorator captures the function being decorated, and wrapper does the actual work using the call_timestamps list from the enclosing scope. Each decorated function gets its own independent rate limit state because each call to rate_limit creates a new call_timestamps list.

Frequently Asked Questions

What exactly makes a function a closure?

A function becomes a closure when it is defined inside another function and references variables from the enclosing function's scope. The closure "closes over" those variables, keeping them alive even after the enclosing function returns. You can check if a function is a closure by inspecting its __closure__ attribute -- if it is not None, the function is a closure.

Why do closures in loops all share the same variable?

This is the most common closure pitfall. When you create closures inside a loop, they all reference the same loop variable, not a copy of it. By the time the closures run, the loop variable has its final value. Fix this by using a default argument: lambda i=i: i captures the current value of i at each iteration.

Are closures slower than regular functions?

The overhead is negligible. Accessing a free variable through a cell object adds one extra pointer dereference compared to accessing a local variable. In practice, this difference is unmeasurable. Closures are used internally by Python for decorators, generators, and comprehensions, so they are highly optimized.

Do closures cause memory leaks?

Closures keep their free variables alive as long as the closure exists, which can prevent garbage collection of those variables. This is rarely a problem in practice, but if your closure captures a large object (like a database connection or a huge list), be aware that the object will not be garbage collected until the closure itself is collected.

Are decorators just closures?

Most decorators are implemented as closures, yes. The decorator function takes the original function as an argument and returns a wrapper closure that adds behavior before or after calling the original. However, decorators can also be implemented as classes with a __call__ method -- the decorator pattern is about the wrapping behavior, not the specific implementation technique.

Conclusion

You have learned how Python closures work from the ground up -- starting with nested functions and the LEGB scoping rule, through the nonlocal keyword for modifying enclosed state, to practical patterns like factory functions, memoization caches, and rate limiters. Closures give you stateful functions without the ceremony of defining a class.

Try refactoring one of your existing classes that has a single method and minimal state into a closure-based factory function. You will be surprised how much simpler the code becomes. For advanced closure patterns, explore Python's functools module, which provides closure-based utilities like lru_cache, partial, and wraps.

Related Articles

• How To Use Python Decorators To Add Logging and Timing to Functions
• How To Use pyproject.toml for Modern Python Project Setup
• How To Use Python coverage.py for Measuring Test Coverage
• How To Use Python pre-commit for Automated Code Quality Hooks
• How To Use Python poetry for Dependency and Project Management
• How To Use Python unittest.mock for Mocking Dependencies
• How To Deploy a Python App to AWS Lambda