From 1deb02959ef822c110c84e944427af2f514b9381 Mon Sep 17 00:00:00 2001
From: Real Python Bot <42617967+realpython-bot@users.noreply.github.com>
Date: Mon, 14 Sep 2026 09:56:45 +0000
Subject: [PATCH] Add materials for the Python statistics fundamentals tutorial

Python Statistics Fundamentals: How to Describe Your Data shipped without
a companion folder. Its 2026 maintenance pass re-ran the tutorial on the
current stack, so add the code here in runnable form.

The tutorial teaches in the REPL. Each of its sections becomes a script
that runs end to end, with the example data shared in datasets.py so a
reader can swap in their own numbers in one place.

Verified on the tutorial's pinned stack (Python 3.14.6, NumPy 2.5.3,
SciPy 1.18.1, pandas 3.0.5, Matplotlib 3.11.2): all 13 scripts run, and
the values match the tutorial, including the NumPy 2 scalar reprs and the
linregress result.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 python-statistics/README.md               |  54 ++++++++
 python-statistics/central_tendency.py     | 154 ++++++++++++++++++++++
 python-statistics/constraints.txt         |   4 +
 python-statistics/correlation.py          |  71 ++++++++++
 python-statistics/dataframes.py           |  64 +++++++++
 python-statistics/datasets.py             |  68 ++++++++++
 python-statistics/plot_bar.py             |  27 ++++
 python-statistics/plot_box.py             |  33 +++++
 python-statistics/plot_heatmap.py         |  36 +++++
 python-statistics/plot_histogram.py       |  30 +++++
 python-statistics/plot_pie.py             |  20 +++
 python-statistics/plot_xy.py              |  31 +++++
 python-statistics/requirements.txt        |   4 +
 python-statistics/summary_statistics.py   |  47 +++++++
 python-statistics/two_dimensional_data.py |  65 +++++++++
 python-statistics/variability.py          | 146 ++++++++++++++++++++
 16 files changed, 854 insertions(+)
 create mode 100644 python-statistics/README.md
 create mode 100644 python-statistics/central_tendency.py
 create mode 100644 python-statistics/constraints.txt
 create mode 100644 python-statistics/correlation.py
 create mode 100644 python-statistics/dataframes.py
 create mode 100644 python-statistics/datasets.py
 create mode 100644 python-statistics/plot_bar.py
 create mode 100644 python-statistics/plot_box.py
 create mode 100644 python-statistics/plot_heatmap.py
 create mode 100644 python-statistics/plot_histogram.py
 create mode 100644 python-statistics/plot_pie.py
 create mode 100644 python-statistics/plot_xy.py
 create mode 100644 python-statistics/requirements.txt
 create mode 100644 python-statistics/summary_statistics.py
 create mode 100644 python-statistics/two_dimensional_data.py
 create mode 100644 python-statistics/variability.py

diff --git a/python-statistics/README.md b/python-statistics/README.md
new file mode 100644
index 0000000000..143a724ca1
--- /dev/null
+++ b/python-statistics/README.md
@@ -0,0 +1,54 @@
+# Python Statistics Fundamentals: How to Describe Your Data
+
+This folder contains supplementary code for the Real Python tutorial on [descriptive statistics in Python](https://realpython.com/python-statistics/). You can run the examples as they are, or continue your learning by experimenting with them on your own data.
+
+The tutorial works through the examples in the REPL. Here, each section of the tutorial is collected into a script you can run end to end, so you can see every result at once and then change the numbers to see what happens.
+
+## Setup
+
+Create and activate a virtual environment, then install the requirements:
+
+```bash
+$ python -m venv venv
+$ source venv/bin/activate
+(venv) $ python -m pip install -r requirements.txt -c constraints.txt
+```
+
+The tutorial's examples were verified against the versions pinned in `constraints.txt`. Newer releases will usually work too, but NumPy 2 changed how scalars are displayed, so output from older versions won't match the tutorial exactly.
+
+## Usage
+
+After creating and activating your virtual environment, and installing the dependencies, you should be able to run each individual file normally:
+
+```bash
+(venv) $ python central_tendency.py
+```
+
+The statistics scripts print their results to the terminal. The plotting scripts open a Matplotlib window instead, one per figure — close each window to move on to the next.
+
+| File | Tutorial section |
+| --- | --- |
+| `datasets.py` | The example data every other script imports |
+| `central_tendency.py` | Measures of Central Tendency |
+| `variability.py` | Measures of Variability |
+| `summary_statistics.py` | Summary of Descriptive Statistics |
+| `correlation.py` | Measures of Correlation Between Pairs of Data |
+| `two_dimensional_data.py` | Working With 2D Data: Axes |
+| `dataframes.py` | Working With 2D Data: DataFrames |
+| `plot_box.py` | Visualizing Data: Box Plots |
+| `plot_histogram.py` | Visualizing Data: Histograms |
+| `plot_pie.py` | Visualizing Data: Pie Charts |
+| `plot_bar.py` | Visualizing Data: Bar Charts |
+| `plot_xy.py` | Visualizing Data: X-Y Plots |
+| `plot_heatmap.py` | Visualizing Data: Heatmaps |
+
+Each script is organized into functions named after the statistic they calculate, and running a file calls them in the order the tutorial introduces them. To work with just one measure, import it instead:
+
+```pycon
+>>> from central_tendency import harmonic_mean
+>>> harmonic_mean()
+```
+
+The example data lives in `datasets.py`, so you can swap in your own numbers in one place and re-run any script against them.
+
+You can find more information and context on the code in [Python Statistics Fundamentals: How to Describe Your Data](https://realpython.com/python-statistics/).
diff --git a/python-statistics/central_tendency.py b/python-statistics/central_tendency.py
new file mode 100644
index 0000000000..9c4b2d1a28
--- /dev/null
+++ b/python-statistics/central_tendency.py
@@ -0,0 +1,154 @@
+"""Measures of central tendency: mean, median, and mode.
+
+Covers the "Measures of Central Tendency" section of
+https://realpython.com/python-statistics/
+"""
+
+import math
+import statistics
+
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from datasets import one_dimensional
+
+
+def mean():
+    """Calculate the arithmetic mean with pure Python, NumPy, and pandas."""
+    x, x_with_nan, y, y_with_nan, z, z_with_nan = one_dimensional()
+
+    print("Pure Python:")
+    print(f"  sum(x) / len(x)          = {sum(x) / len(x)}")
+    print(f"  statistics.mean(x)       = {statistics.mean(x)}")
+    print(f"  statistics.fmean(x)      = {statistics.fmean(x)}")
+
+    # A single nan in the data makes the result nan.
+    print(f"  statistics.mean(x_with_nan)  = {statistics.mean(x_with_nan)}")
+
+    print("NumPy:")
+    print(f"  np.mean(y)               = {np.mean(y)}")
+    print(f"  y.mean()                 = {y.mean()}")
+    print(f"  np.mean(y_with_nan)      = {np.mean(y_with_nan)}")
+    # Use the nan-safe variant to ignore missing values instead.
+    print(f"  np.nanmean(y_with_nan)   = {np.nanmean(y_with_nan)}")
+
+    print("pandas:")
+    # pandas skips nan values by default.
+    print(f"  z.mean()                 = {z.mean()}")
+    print(f"  z_with_nan.mean()        = {z_with_nan.mean()}")
+
+
+def weighted_mean():
+    """Calculate the weighted mean, where each item has its own weight."""
+    x = [8.0, 1, 2.5, 4, 28.0]
+    w = [0.1, 0.2, 0.3, 0.25, 0.15]
+
+    wmean = sum(w_ * x_ for (x_, w_) in zip(x, w)) / sum(w)
+    print(f"Pure Python: {wmean}")
+
+    y, z, w_array = np.array(x), pd.Series(x), np.array(w)
+    print(f"np.average(y, weights=w) = {np.average(y, weights=w_array)}")
+    print(f"np.average(z, weights=w) = {np.average(z, weights=w_array)}")
+    print(f"(w * y).sum() / w.sum()  = {(w_array * y).sum() / w_array.sum()}")
+
+
+def harmonic_mean():
+    """Calculate the harmonic mean, the reciprocal of the mean reciprocal."""
+    x, _, y, _, z, _ = one_dimensional()
+
+    hmean = len(x) / sum(1 / item for item in x)
+    print(f"Pure Python:                 {hmean}")
+    print(f"statistics.harmonic_mean(x): {statistics.harmonic_mean(x)}")
+    print(f"scipy.stats.hmean(y):        {scipy.stats.hmean(y)}")
+    print(f"scipy.stats.hmean(z):        {scipy.stats.hmean(z)}")
+
+    # A negative value has no harmonic mean.
+    try:
+        statistics.harmonic_mean([1, 2, -2])
+    except statistics.StatisticsError as error:
+        print(f"harmonic_mean([1, 2, -2]) raises StatisticsError: {error}")
+
+
+def geometric_mean():
+    """Calculate the geometric mean, the nth root of the product."""
+    x, _, y, _, z, _ = one_dimensional()
+
+    gmean = 1
+    for item in x:
+        gmean *= item
+    gmean **= 1 / len(x)
+
+    print(f"Pure Python:                  {gmean}")
+    print(f"statistics.geometric_mean(x): {statistics.geometric_mean(x)}")
+    print(f"scipy.stats.gmean(y):         {scipy.stats.gmean(y)}")
+    print(f"scipy.stats.gmean(z):         {scipy.stats.gmean(z)}")
+
+
+def median():
+    """Find the median, the middle value of the sorted data."""
+    x, x_with_nan, y, y_with_nan, z, z_with_nan = one_dimensional()
+
+    n = len(x)
+    if n % 2:
+        median_ = sorted(x)[round(0.5 * (n - 1))]
+    else:
+        x_ord, index = sorted(x), round(0.5 * n)
+        median_ = 0.5 * (x_ord[index - 1] + x_ord[index])
+    print(f"Pure Python:          {median_}")
+
+    print(f"statistics.median(x): {statistics.median(x)}")
+    # With an even number of items, the median averages the two middle values.
+    print(f"median(x[:-1]):       {statistics.median(x[:-1])}")
+    print(f"median_low(x[:-1]):   {statistics.median_low(x[:-1])}")
+    print(f"median_high(x[:-1]):  {statistics.median_high(x[:-1])}")
+
+    print(f"np.median(y):         {np.median(y)}")
+    print(f"np.nanmedian(y_with_nan): {np.nanmedian(y_with_nan)}")
+    print(f"z.median():           {z.median()}")
+    print(f"z_with_nan.median():  {z_with_nan.median()}")
+
+
+def mode():
+    """Find the mode, the value that appears most often."""
+    u = [2, 3, 2, 8, 12]
+    v = [12, 15, 12, 15, 21, 15, 12]
+
+    mode_ = max((u.count(item), item) for item in set(u))[1]
+    print(f"Pure Python:            {mode_}")
+    print(f"statistics.mode(u):     {statistics.mode(u)}")
+    print(f"statistics.multimode(u): {statistics.multimode(u)}")
+
+    # multimode() returns every modal value, mode() only the first.
+    print(f"statistics.mode(v):     {statistics.mode(v)}")
+    print(f"statistics.multimode(v): {statistics.multimode(v)}")
+
+    u_array, v_array = np.array(u), np.array(v)
+    print(f"scipy.stats.mode(u):    {scipy.stats.mode(u_array)}")
+
+    result = scipy.stats.mode(v_array)
+    print(f"scipy.stats.mode(v):    {result}")
+    print(f"  .mode  = {result.mode}")
+    print(f"  .count = {result.count}")
+
+    # pandas returns a Series, so it can report several modal values at once.
+    u_series = pd.Series(u)
+    v_series = pd.Series(v)
+    w_series = pd.Series([2, 2, math.nan])
+    print(f"u.mode():\n{u_series.mode()}")
+    print(f"v.mode():\n{v_series.mode()}")
+    print(f"w.mode():\n{w_series.mode()}")
+
+
+if __name__ == "__main__":
+    for section in (
+        mean,
+        weighted_mean,
+        harmonic_mean,
+        geometric_mean,
+        median,
+        mode,
+    ):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()
diff --git a/python-statistics/constraints.txt b/python-statistics/constraints.txt
new file mode 100644
index 0000000000..9a278490b4
--- /dev/null
+++ b/python-statistics/constraints.txt
@@ -0,0 +1,4 @@
+numpy==2.5.3
+scipy==1.18.1
+pandas==3.0.5
+matplotlib==3.11.2
diff --git a/python-statistics/correlation.py b/python-statistics/correlation.py
new file mode 100644
index 0000000000..2aad5267f0
--- /dev/null
+++ b/python-statistics/correlation.py
@@ -0,0 +1,71 @@
+"""Measures of correlation between pairs of data.
+
+Covers the "Measures of Correlation Between Pairs of Data" section of
+https://realpython.com/python-statistics/
+"""
+
+import numpy as np
+import scipy.stats
+
+from datasets import correlation_data
+
+
+def covariance():
+    """Calculate the covariance, which measures joint variability."""
+    x, y, x_, y_, x__, y__ = correlation_data()
+
+    n = len(x)
+    mean_x, mean_y = sum(x) / n, sum(y) / n
+    cov_xy = sum((x[k] - mean_x) * (y[k] - mean_y) for k in range(n)) / (n - 1)
+    print(f"Pure Python:            {cov_xy}")
+
+    # np.cov() returns the covariance matrix, not a single number. Its
+    # diagonal holds the variances and its off-diagonal the covariance.
+    cov_matrix = np.cov(x_, y_)
+    print(f"np.cov(x_, y_) =\n{cov_matrix}")
+    print(f"  variance of x: {x_.var(ddof=1)}")
+    print(f"  variance of y: {y_.var(ddof=1)}")
+    print(f"  covariance:    {cov_matrix[0, 1]}")
+
+    print(f"x__.cov(y__):           {x__.cov(y__)}")
+
+
+def correlation_coefficient():
+    """Calculate Pearson's r, the normalized measure of correlation."""
+    x, y, x_, y_, x__, y__ = correlation_data()
+
+    n = len(x)
+    mean_x, mean_y = sum(x) / n, sum(y) / n
+    var_x = sum((item - mean_x) ** 2 for item in x) / (n - 1)
+    var_y = sum((item - mean_y) ** 2 for item in y) / (n - 1)
+    std_x, std_y = var_x**0.5, var_y**0.5
+    cov_xy = sum((x[k] - mean_x) * (y[k] - mean_y) for k in range(n)) / (n - 1)
+    print(f"Pure Python:            {cov_xy / (std_x * std_y)}")
+
+    # pearsonr() returns the coefficient and the p-value together.
+    r, p = scipy.stats.pearsonr(x_, y_)
+    print(f"scipy.stats.pearsonr(): r={r}, p={p}")
+
+    corr_matrix = np.corrcoef(x_, y_)
+    print(f"np.corrcoef(x_, y_) =\n{corr_matrix}")
+    print(f"  r = {corr_matrix[0, 1]}")
+
+    print(f"x__.corr(y__):          {x__.corr(y__)}")
+
+
+def linear_regression():
+    """Fit a regression line, which also reports the correlation."""
+    *_, x_, y_, __, ___ = correlation_data()
+
+    result = scipy.stats.linregress(x_, y_)
+    print(result)
+    print(f"  slope     = {result.slope}")
+    print(f"  intercept = {result.intercept}")
+    print(f"  rvalue    = {result.rvalue}")
+
+
+if __name__ == "__main__":
+    for section in (covariance, correlation_coefficient, linear_regression):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()
diff --git a/python-statistics/dataframes.py b/python-statistics/dataframes.py
new file mode 100644
index 0000000000..36cdd7670c
--- /dev/null
+++ b/python-statistics/dataframes.py
@@ -0,0 +1,64 @@
+"""Statistics with pandas DataFrame objects.
+
+Covers the "DataFrames" section of https://realpython.com/python-statistics/
+"""
+
+import pandas as pd
+
+from datasets import two_dimensional
+
+
+def build_dataframe():
+    """Return the labeled DataFrame used throughout this script."""
+    row_names = ["first", "second", "third", "fourth", "fifth"]
+    col_names = ["A", "B", "C"]
+    return pd.DataFrame(two_dimensional(), index=row_names, columns=col_names)
+
+
+def by_column_and_row():
+    """Calculate statistics across columns and across rows."""
+    df = build_dataframe()
+    print(f"df =\n{df}\n")
+
+    # DataFrame methods default to working column by column.
+    print(f"df.mean():\n{df.mean()}\n")
+    print(f"df.var():\n{df.var()}\n")
+
+    # Pass axis=1 to work row by row instead.
+    print(f"df.mean(axis=1):\n{df.mean(axis=1)}\n")
+    print(f"df.var(axis=1):\n{df.var(axis=1)}")
+
+
+def single_column():
+    """Pull one column out as a Series and work with it directly."""
+    df = build_dataframe()
+
+    print(f"df['A']:\n{df['A']}\n")
+    print(f"df['A'].mean(): {df['A'].mean()}")
+    print(f"df['A'].var():  {df['A'].var()}")
+
+
+def to_numpy():
+    """Drop the labels and get back a plain NumPy array."""
+    df = build_dataframe()
+
+    print(f"df.values:\n{df.values}\n")
+    print(f"df.to_numpy():\n{df.to_numpy()}")
+
+
+def describe():
+    """Summarize every column at once."""
+    df = build_dataframe()
+
+    print(f"df.describe():\n{df.describe()}\n")
+
+    # The summary is a DataFrame too, so .at[] reaches a single value.
+    print(f"df.describe().at['mean', 'A']: {df.describe().at['mean', 'A']}")
+    print(f"df.describe().at['50%', 'B']:  {df.describe().at['50%', 'B']}")
+
+
+if __name__ == "__main__":
+    for section in (by_column_and_row, single_column, to_numpy, describe):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()
diff --git a/python-statistics/datasets.py b/python-statistics/datasets.py
new file mode 100644
index 0000000000..a360725ea5
--- /dev/null
+++ b/python-statistics/datasets.py
@@ -0,0 +1,68 @@
+"""Example datasets used throughout the tutorial.
+
+Each function returns the data for one section of
+https://realpython.com/python-statistics/ so that the other scripts in this
+folder can stay short and focused on the statistics themselves.
+"""
+
+import math
+
+import numpy as np
+import pandas as pd
+
+
+def one_dimensional():
+    """Return the small 1D dataset used for central tendency and variability.
+
+    Returns a plain list, a NumPy array, and a pandas Series, each in a
+    version with and without a ``nan`` value.
+    """
+    x = [8.0, 1, 2.5, 4, 28.0]
+    x_with_nan = [8.0, 1, 2.5, math.nan, 4, 28.0]
+
+    y, y_with_nan = np.array(x), np.array(x_with_nan)
+    z, z_with_nan = pd.Series(x), pd.Series(x_with_nan)
+
+    return x, x_with_nan, y, y_with_nan, z, z_with_nan
+
+
+def percentile_data():
+    """Return the nine-value dataset used in the percentiles section."""
+    x = [-5.0, -1.1, 0.1, 2.0, 8.0, 12.8, 21.0, 25.8, 41.0]
+    y = np.array(x)
+    y_with_nan = np.insert(y, 2, np.nan)
+
+    return x, y, y_with_nan
+
+
+def correlation_data():
+    """Return the paired dataset used for covariance and correlation."""
+    x = list(range(-10, 11))
+    y = [0, 2, 2, 2, 2, 3, 3, 6, 7, 4, 7, 6, 6, 9, 4, 5, 5, 10, 11, 12, 14]
+
+    x_, y_ = np.array(x), np.array(y)
+    x__, y__ = pd.Series(x_), pd.Series(y_)
+
+    return x, y, x_, y_, x__, y__
+
+
+def two_dimensional():
+    """Return the 5x3 array used in the sections on axes and DataFrames."""
+    return np.array(
+        [
+            [1, 1, 1],
+            [2, 3, 1],
+            [4, 9, 2],
+            [8, 27, 4],
+            [16, 1, 1],
+        ]
+    )
+
+
+if __name__ == "__main__":
+    x, x_with_nan, y, y_with_nan, z, z_with_nan = one_dimensional()
+    print(f"{x = }")
+    print(f"{x_with_nan = }")
+    print(f"y =\n{y}\n")
+    print(f"z =\n{z}\n")
+    print(f"2D array a =\n{two_dimensional()}")
diff --git a/python-statistics/plot_bar.py b/python-statistics/plot_bar.py
new file mode 100644
index 0000000000..5d3937eb0e
--- /dev/null
+++ b/python-statistics/plot_bar.py
@@ -0,0 +1,27 @@
+"""Draw a bar chart with error bars.
+
+Covers the "Bar Charts" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+
+def main():
+    np.random.seed(seed=0)
+    x = np.arange(21)
+    y = np.random.randint(21, size=21)
+    # Error bar sizes must be non-negative, so take the absolute values.
+    err = np.abs(np.random.randn(21))
+
+    fig, ax = plt.subplots()
+    ax.bar(x, y, yerr=err)
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/plot_box.py b/python-statistics/plot_box.py
new file mode 100644
index 0000000000..134088bfd3
--- /dev/null
+++ b/python-statistics/plot_box.py
@@ -0,0 +1,33 @@
+"""Draw a box plot of three random datasets.
+
+Covers the "Box Plots" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+
+def main():
+    np.random.seed(seed=0)
+    x = np.random.randn(1000)
+    y = np.random.randn(100)
+    z = np.random.randn(10)
+
+    fig, ax = plt.subplots()
+    ax.boxplot(
+        (x, y, z),
+        orientation="horizontal",
+        showmeans=True,
+        meanline=True,
+        tick_labels=("x", "y", "z"),
+        patch_artist=True,
+        medianprops={"linewidth": 2, "color": "purple"},
+        meanprops={"linewidth": 2, "color": "red"},
+    )
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/plot_heatmap.py b/python-statistics/plot_heatmap.py
new file mode 100644
index 0000000000..45f615a6c7
--- /dev/null
+++ b/python-statistics/plot_heatmap.py
@@ -0,0 +1,36 @@
+"""Draw heatmaps of the covariance and correlation matrices.
+
+Covers the "Heatmaps" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+
+def draw_heatmap(matrix):
+    """Draw one 2x2 matrix as a heatmap with its values written on top."""
+    fig, ax = plt.subplots()
+    ax.imshow(matrix)
+    ax.grid(False)
+    ax.xaxis.set(ticks=(0, 1), ticklabels=("x", "y"))
+    ax.yaxis.set(ticks=(0, 1), ticklabels=("x", "y"))
+    ax.set_ylim(1.5, -0.5)
+    for i in range(2):
+        for j in range(2):
+            ax.text(j, i, matrix[i, j], ha="center", va="center", color="w")
+    plt.show()
+
+
+def main():
+    np.random.seed(seed=0)
+    x = np.arange(21)
+    y = 5 + 2 * x + 2 * np.random.randn(21)
+
+    draw_heatmap(np.cov(x, y).round(decimals=2))
+    draw_heatmap(np.corrcoef(x, y).round(decimals=2))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/plot_histogram.py b/python-statistics/plot_histogram.py
new file mode 100644
index 0000000000..01179f9fe9
--- /dev/null
+++ b/python-statistics/plot_histogram.py
@@ -0,0 +1,30 @@
+"""Draw ordinary and cumulative histograms.
+
+Covers the "Histograms" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+
+def main():
+    np.random.seed(seed=0)
+    x = np.random.randn(1000)
+
+    # np.histogram() gives you the same numbers that .hist() draws.
+    hist, bin_edges = np.histogram(x, bins=10)
+    print(f"{hist = }")
+    print(f"{bin_edges = }")
+
+    for cumulative in (False, True):
+        fig, ax = plt.subplots()
+        ax.hist(x, bin_edges, cumulative=cumulative)
+        ax.set_xlabel("x")
+        ax.set_ylabel("Frequency")
+        plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/plot_pie.py b/python-statistics/plot_pie.py
new file mode 100644
index 0000000000..8b42b61376
--- /dev/null
+++ b/python-statistics/plot_pie.py
@@ -0,0 +1,20 @@
+"""Draw a pie chart of three values.
+
+Covers the "Pie Charts" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+
+plt.style.use("ggplot")
+
+
+def main():
+    x, y, z = 128, 256, 1024
+
+    fig, ax = plt.subplots()
+    ax.pie((x, y, z), labels=("x", "y", "z"), autopct="%1.1f%%")
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/plot_xy.py b/python-statistics/plot_xy.py
new file mode 100644
index 0000000000..9723710abb
--- /dev/null
+++ b/python-statistics/plot_xy.py
@@ -0,0 +1,31 @@
+"""Draw an x-y plot with its regression line.
+
+Covers the "X-Y Plots" section of https://realpython.com/python-statistics/
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.stats
+
+plt.style.use("ggplot")
+
+
+def main():
+    np.random.seed(seed=0)
+    x = np.arange(21)
+    y = 5 + 2 * x + 2 * np.random.randn(21)
+
+    slope, intercept, r, *__ = scipy.stats.linregress(x, y)
+    line = f"Regression line: y={intercept:.2f}+{slope:.2f}x, r={r:.2f}"
+
+    fig, ax = plt.subplots()
+    ax.plot(x, y, linewidth=0, marker="s", label="Data points")
+    ax.plot(x, intercept + slope * x, label=line)
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    ax.legend(facecolor="white")
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python-statistics/requirements.txt b/python-statistics/requirements.txt
new file mode 100644
index 0000000000..7bd7c445cc
--- /dev/null
+++ b/python-statistics/requirements.txt
@@ -0,0 +1,4 @@
+numpy
+scipy
+pandas
+matplotlib
diff --git a/python-statistics/summary_statistics.py b/python-statistics/summary_statistics.py
new file mode 100644
index 0000000000..cb3e4a0da6
--- /dev/null
+++ b/python-statistics/summary_statistics.py
@@ -0,0 +1,47 @@
+"""Summarize a whole dataset in a single call.
+
+Covers the "Summary of Descriptive Statistics" section of
+https://realpython.com/python-statistics/
+"""
+
+import pandas as pd
+import scipy.stats
+
+from datasets import percentile_data
+
+
+def scipy_describe():
+    """Summarize an array with scipy.stats.describe()."""
+    _, y, _ = percentile_data()
+
+    result = scipy.stats.describe(y, ddof=1, bias=False)
+    print(result)
+
+    # Every field is also available as an attribute.
+    print(f"  nobs      = {result.nobs}")
+    print(f"  min       = {result.minmax[0]}")
+    print(f"  max       = {result.minmax[1]}")
+    print(f"  mean      = {result.mean}")
+    print(f"  variance  = {result.variance}")
+    print(f"  skewness  = {result.skewness}")
+    print(f"  kurtosis  = {result.kurtosis}")
+
+
+def pandas_describe():
+    """Summarize a Series with .describe()."""
+    _, y, _ = percentile_data()
+    z = pd.Series(y)
+
+    result = z.describe()
+    print(result)
+
+    # The result is itself a Series, so you index it by label.
+    for label in ("mean", "std", "min", "max", "25%", "50%", "75%"):
+        print(f"  {label:<4} = {result[label]}")
+
+
+if __name__ == "__main__":
+    for section in (scipy_describe, pandas_describe):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()
diff --git a/python-statistics/two_dimensional_data.py b/python-statistics/two_dimensional_data.py
new file mode 100644
index 0000000000..c6c29caa8d
--- /dev/null
+++ b/python-statistics/two_dimensional_data.py
@@ -0,0 +1,65 @@
+"""Statistics over 2D data, and what the axis argument does.
+
+Covers the "Working With 2D Data" section of
+https://realpython.com/python-statistics/
+"""
+
+import numpy as np
+import scipy.stats
+
+from datasets import two_dimensional
+
+
+def whole_array():
+    """Apply statistics to every value in the array at once."""
+    a = two_dimensional()
+    print(f"a =\n{a}")
+
+    print(f"np.mean(a):             {np.mean(a)}")
+    print(f"a.mean():               {a.mean()}")
+    print(f"np.median(a):           {np.median(a)}")
+    print(f"a.var(ddof=1):          {a.var(ddof=1)}")
+
+
+def along_axes():
+    """Apply statistics column by column (axis=0) or row by row (axis=1)."""
+    a = two_dimensional()
+
+    # axis=0 collapses the rows, leaving one result per column.
+    print(f"np.mean(a, axis=0):     {np.mean(a, axis=0)}")
+    print(f"np.median(a, axis=0):   {np.median(a, axis=0)}")
+    print(f"a.var(axis=0, ddof=1):  {a.var(axis=0, ddof=1)}")
+
+    # axis=1 collapses the columns, leaving one result per row.
+    print(f"np.mean(a, axis=1):     {np.mean(a, axis=1)}")
+    print(f"np.median(a, axis=1):   {np.median(a, axis=1)}")
+    print(f"a.var(axis=1, ddof=1):  {a.var(axis=1, ddof=1)}")
+
+    # SciPy functions take the same argument, and default to axis=0.
+    print(f"scipy.stats.gmean(a):   {scipy.stats.gmean(a)}")
+    print(f"gmean(a, axis=1):       {scipy.stats.gmean(a, axis=1)}")
+    # axis=None flattens the array first.
+    print(f"gmean(a, axis=None):    {scipy.stats.gmean(a, axis=None)}")
+
+
+def describe_by_axis():
+    """Summarize a 2D array, flattened and along each axis."""
+    a = two_dimensional()
+
+    print("axis=None (the whole array):")
+    print(f"  {scipy.stats.describe(a, axis=None, ddof=1, bias=False)}")
+
+    print("axis=0 (per column, the default):")
+    print(f"  {scipy.stats.describe(a, ddof=1, bias=False)}")
+
+    print("axis=1 (per row):")
+    result = scipy.stats.describe(a, axis=1, ddof=1, bias=False)
+    print(f"  {result}")
+    print(f"  .mean = {result.mean}")
+
+
+if __name__ == "__main__":
+    for section in (whole_array, along_axes, describe_by_axis):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()
diff --git a/python-statistics/variability.py b/python-statistics/variability.py
new file mode 100644
index 0000000000..d286fb2c0d
--- /dev/null
+++ b/python-statistics/variability.py
@@ -0,0 +1,146 @@
+"""Measures of variability: variance, standard deviation, skewness, and more.
+
+Covers the "Measures of Variability" section of
+https://realpython.com/python-statistics/
+"""
+
+import statistics
+
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from datasets import one_dimensional, percentile_data
+
+
+def variance():
+    """Calculate the sample variance, the mean squared deviation."""
+    x, x_with_nan, y, y_with_nan, z, z_with_nan = one_dimensional()
+
+    n = len(x)
+    mean_ = sum(x) / n
+    var_ = sum((item - mean_) ** 2 for item in x) / (n - 1)
+    print(f"Pure Python:            {var_}")
+    print(f"statistics.variance(x): {statistics.variance(x)}")
+
+    # ddof=1 asks NumPy for the sample variance rather than the population one.
+    print(f"np.var(y, ddof=1):      {np.var(y, ddof=1)}")
+    print(f"y.var(ddof=1):          {y.var(ddof=1)}")
+    print(f"np.nanvar(y_with_nan, ddof=1): {np.nanvar(y_with_nan, ddof=1)}")
+
+    # pandas defaults to ddof=1 and skips nan values.
+    print(f"z.var(ddof=1):          {z.var(ddof=1)}")
+    print(f"z_with_nan.var(ddof=1): {z_with_nan.var(ddof=1)}")
+
+
+def standard_deviation():
+    """Calculate the standard deviation, the square root of the variance."""
+    x, _, y, y_with_nan, z, z_with_nan = one_dimensional()
+
+    var_ = statistics.variance(x)
+    print(f"var_ ** 0.5:            {var_**0.5}")
+    print(f"statistics.stdev(x):    {statistics.stdev(x)}")
+
+    print(f"np.std(y, ddof=1):      {np.std(y, ddof=1)}")
+    print(f"y.std(ddof=1):          {y.std(ddof=1)}")
+    print(f"np.nanstd(y_with_nan, ddof=1): {np.nanstd(y_with_nan, ddof=1)}")
+
+    print(f"z.std(ddof=1):          {z.std(ddof=1)}")
+    print(f"z_with_nan.std(ddof=1): {z_with_nan.std(ddof=1)}")
+
+
+def skewness():
+    """Measure skewness, how asymmetric the data is around its mean."""
+    x, x_with_nan, y, y_with_nan, z, z_with_nan = one_dimensional()
+
+    n = len(x)
+    mean_ = sum(x) / n
+    var_ = sum((item - mean_) ** 2 for item in x) / (n - 1)
+    std_ = var_**0.5
+    skew_ = (
+        sum((item - mean_) ** 3 for item in x)
+        * n
+        / ((n - 1) * (n - 2) * std_**3)
+    )
+    print(f"Pure Python:            {skew_}")
+
+    # bias=False applies the correction for statistical bias.
+    print(
+        f"scipy.stats.skew(y, bias=False): {scipy.stats.skew(y, bias=False)}"
+    )
+    print(f"z.skew():               {z.skew()}")
+    print(f"z_with_nan.skew():      {z_with_nan.skew()}")
+
+
+def percentiles():
+    """Find percentiles and quantiles, the cut points that divide the data."""
+    x, y, y_with_nan = percentile_data()
+
+    print(f"quantiles(x, n=2):      {statistics.quantiles(x, n=2)}")
+    print(
+        "quantiles(x, n=4, inclusive): "
+        f"{statistics.quantiles(x, n=4, method='inclusive')}"
+    )
+
+    print(f"np.percentile(y, 5):    {np.percentile(y, 5)}")
+    print(f"np.percentile(y, 95):   {np.percentile(y, 95)}")
+    print(f"np.percentile(y, [25, 50, 75]): {np.percentile(y, [25, 50, 75])}")
+    print(
+        "np.nanpercentile(y_with_nan, [25, 50, 75]): "
+        f"{np.nanpercentile(y_with_nan, [25, 50, 75])}"
+    )
+
+    # quantile() takes fractions between 0 and 1 instead of percentages.
+    print(f"np.quantile(y, 0.05):   {np.quantile(y, 0.05)}")
+    print(
+        f"np.quantile(y, [0.25, 0.5, 0.75]): {np.quantile(y, [0.25, 0.5, 0.75])}"
+    )
+
+    z, z_with_nan = pd.Series(y), pd.Series(y_with_nan)
+    print(f"z.quantile(0.05):       {z.quantile(0.05)}")
+    print(f"z.quantile([0.25, 0.5, 0.75]):\n{z.quantile([0.25, 0.5, 0.75])}")
+    print(
+        "z_with_nan.quantile([0.25, 0.5, 0.75]):\n"
+        f"{z_with_nan.quantile([0.25, 0.5, 0.75])}"
+    )
+
+
+def ranges():
+    """Calculate the range and the interquartile range."""
+    x, y, y_with_nan = percentile_data()
+    z, z_with_nan = pd.Series(y), pd.Series(y_with_nan)
+
+    # np.ptp() returns nan when the data contains nan, whether you pass a
+    # NumPy array or a pandas Series.
+    print(f"np.ptp(y):              {np.ptp(y)}")
+    print(f"np.ptp(z):              {np.ptp(z)}")
+    print(f"np.ptp(y_with_nan):     {np.ptp(y_with_nan)}")
+    print(f"np.ptp(z_with_nan):     {np.ptp(z_with_nan)}")
+
+    print(f"np.amax(y) - np.amin(y): {np.amax(y) - np.amin(y)}")
+    print(
+        "np.nanmax(y_with_nan) - np.nanmin(y_with_nan): "
+        f"{np.nanmax(y_with_nan) - np.nanmin(y_with_nan)}"
+    )
+    print(f"z.max() - z.min():      {z.max() - z.min()}")
+
+    # The interquartile range is the distance between the first and third
+    # quartiles, which ignores outliers at both ends.
+    quartiles = np.quantile(y, [0.25, 0.75])
+    print(f"Interquartile range (NumPy):  {quartiles[1] - quartiles[0]}")
+
+    quartiles = z.quantile([0.25, 0.75])
+    print(f"Interquartile range (pandas): {quartiles[0.75] - quartiles[0.25]}")
+
+
+if __name__ == "__main__":
+    for section in (
+        variance,
+        standard_deviation,
+        skewness,
+        percentiles,
+        ranges,
+    ):
+        title = section.__name__.replace("_", " ").title()
+        print(f"\n{title}\n{'-' * len(title)}")
+        section()