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>

Author: realpython-bot

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/).

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()

python-statistics/constraints.txt

@@ -0,0 +1,4 @@
+numpy==2.5.3
+scipy==1.18.1
+pandas==3.0.5
+matplotlib==3.11.2

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()

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()