diff --git a/README.md b/README.md
index 0758dd0..dd85a79 100644
--- a/README.md
+++ b/README.md
@@ -76,21 +76,17 @@ my_timeseries = ...
 # of the first sample after equilibration, g is the statistical
 # inefficiency of the equilibrated sample, and ess is the effective sample
 # size of the equilibrated sample.
-idx, g, ess = red.detect_equilibration_window(my_timeseries,
-                                              method="min_sse",
-                                              plot=True)
+idx, g, ess = red.detect_equilibration_window(my_timeseries, method="min_sse", plot=True)
 
 # Alternatively, use Geyer's initial convex sequence method to account
 # for autocorrelation.
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
-                                                method="min_sse",
-                                                plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="min_sse", plot=True)
 
 # We can also determine equilibration in the same way as in
 # pymbar.timeseries.detect_equilibration(my_timeseries, fast=False)
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
-                                                method="max_ess",
-                                                sequence_estimator="positive")
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, method="max_ess", sequence_estimator="positive"
+)
 ```
 
 #### Uncertainty Quantification
@@ -99,7 +95,6 @@ idx, g, ess = red.detect_equilibration_init_seq(my_timeseries,
 # Estimate the 95 % confidence interval, accounting for autocorrelation using Geyer's initial
 # convex sequence method.
 ci_95 = red.get_conf_int_init_seq(my_timeseries, alpha_two_tailed=0.05)
-
 ```
 
 For more examples, see the [documentation](https://fjclark.github.io/red/latest/examples/).
diff --git a/docs/changelog.md b/docs/changelog.md
index b3488fa..cc10d7c 100644
--- a/docs/changelog.md
+++ b/docs/changelog.md
@@ -2,6 +2,7 @@
 
 ## 0.2.0 - 2026-08-06
 
+- Make ruff stricter about docstrings and fix them up.
 - Switch to from conda, make, and mypy to uv, just, and ty. Modernise type hinting and switch to Python > 3.11.
 - Drop polyfill.io due to security issues.
 
diff --git a/docs/examples.md b/docs/examples.md
index 1a9577a..5470f39 100644
--- a/docs/examples.md
+++ b/docs/examples.md
@@ -26,7 +26,9 @@ For all examples, `my_timeseries` should be a numpy array with shape `(n_samples
 To use any of Geyer's initial sequence methods ([Geyer, 1992](https://www.jstor.org/stable/2246094)), you can specify the "sequence_estimator" to be "initial_positive" (the least strict), "initial_monotone", or "initial_convex" (the strictest):
 
 ```python
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, sequence_estimator="initial_convex", plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, sequence_estimator="initial_convex", plot=True
+)
 my_truncated_timeseries = my_timeseries[idx:]
 ```
 To use Chodera's method of simply truncating the autocovariance series at the first negative value ([Chodera, 2016](https://pubs.acs.org/doi/full/10.1021/acs.jctc.5b00784)), you can specify the "sequence estimator" to be "positive".
@@ -36,7 +38,9 @@ To use Chodera's method of simply truncating the autocovariance series at the fi
 When using window methods, you can either specify a fixed window size, or a window size function which computes the window size as a function of the number of data points (which decreases as the truncation point increases). These are specified via `window_size` and `window_size_fn`, respectively (one must be specified and the other must be `None`). The default window size function is `lambda x: round(x**0.5)` - explicitly:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=None, window_size_fn=lambda x: round(x**0.5), plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=None, window_size_fn=lambda x: round(x**0.5), plot=True
+)
 # This is equivalent to:
 idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True)
 ```
@@ -44,7 +48,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True)
 To use a window size of 10:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=10, window_size_fn = None, plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=10, window_size_fn=None, plot=True
+)
 ```
 
 You can also play with the kernel function used in the window method by specifying the `kernel` argument. You should supply the function directly - the default is `np.bartlett`.
@@ -54,7 +60,9 @@ You can also play with the kernel function used in the window method by specifyi
 To use White's original Marginal Standard Error Rule ([White, 1997](https://journals.sagepub.com/doi/abs/10.1177/003754979706900601)), you can use the window method with a window size of 1:
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=1, window_size_fn=None, plot=True)
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, window_size=1, window_size_fn=None, plot=True
+)
 ```
 
 ### Maximum Effective Sample Size and Chodera's Method
@@ -62,7 +70,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, window_size=1, wind
 To select the truncation point according to the maximum effective sample size (instead of the minimum squared standard error), you can specify the `method` argument to be "max_ess". To use Chodera's method ([Chodera, 2016](https://pubs.acs.org/doi/full/10.1021/acs.jctc.5b00784)) as implemented in `pymbar.timeseries`, you can specify the `sequence_estimator` to be "positive":
 
 ```python
-idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="max_ess", sequence_estimator="positive", plot=True)
+idx, g, ess = red.detect_equilibration_init_seq(
+    my_timeseries, method="max_ess", sequence_estimator="positive", plot=True
+)
 # Equivalent to pymbar.timeseries.detect_equilibration(my_timeseries, fast=False)
 ```
 
@@ -71,7 +81,9 @@ idx, g, ess = red.detect_equilibration_init_seq(my_timeseries, method="max_ess",
 To save a plot showing the (block-averaged) time series and variance of the mean/ effective sample size against truncation time, simply specify `plot=True` and, optionally, specify a name for the plot with `plot_name`. This works for either of the equilbration detection functions.
 
 ```python
-idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True, plot_name="my_equilibration_plot.png")
+idx, g, ess = red.detect_equilibration_window(
+    my_timeseries, plot=True, plot_name="my_equilibration_plot.png"
+)
 ```
 
 ## Estimating Uncertainty
@@ -79,7 +91,9 @@ idx, g, ess = red.detect_equilibration_window(my_timeseries, plot=True, plot_nam
 To calculate uncertainty, we recommend using Geyer's initial convex sequence method ([Geyer, 1992](https://www.jstor.org/stable/2246094)), which is the default for [`get_conf_int_init_seq`][red.confidence_intervals.get_conf_int_init_seq]. For example, to estimate a 95 % confidence interval:
 
 ```python
-ci_95 = red.get_conf_int_init_seq(my_timeseries, sequence_estimator="initial_convex", alpha_two_tailed=0.05)
+ci_95 = red.get_conf_int_init_seq(
+    my_timeseries, sequence_estimator="initial_convex", alpha_two_tailed=0.05
+)
 ```
 
 This function has a similar interface to [`detect_equilibration_init_seq`][red.equilibration.detect_equilibration_init_seq], and you can specify the "sequence_estimator" in the same way. Note that this assumes we have a reasonable effective sample size, and hence that the means are approximately normally distributed by the central limit theorem.
diff --git a/pyproject.toml b/pyproject.toml
index 44f80cf..3d18f87 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -81,10 +81,18 @@ line-length = 100
 [tool.ruff.lint]
 ignore = ["PLR", "PLW", "C901"]
 select = ["B","C","E","F","W","B9"]
+# Enforce PEP257 docstring conventions (pydocstyle). The numpy convention matches the
+# numpy-style docstrings used throughout the package; it disables D212/D413/D415/D417 by
+# default, so we re-enable those explicitly. D213 is intentionally left disabled by the
+# numpy convention (D212 and D213 are mutually exclusive).
+extend-select = ["D", "D212", "D413", "D415", "D417"]
+
+[tool.ruff.lint.pydocstyle]
+convention = "numpy"
 
 [tool.ruff.lint.per-file-ignores]
 "__init__.py" = ["F401"]
-"red/tests/*.py" = ["F401", "F811"]
+"red/tests/*.py" = ["F401", "F811", "D"]
 
 [tool.setuptools]
 # This subkey is a beta stage development and keys may change in the future, see https://setuptools.pypa.io/en/latest/userguide/pyproject_config.html for more details
diff --git a/red/__init__.py b/red/__init__.py
index 3852470..6b83bf8 100644
--- a/red/__init__.py
+++ b/red/__init__.py
@@ -1,4 +1,4 @@
-"""Robust Equilibration Detection"""
+"""Robust Equilibration Detection."""
 
 from ._version import __version__
 from .confidence_intervals import get_conf_int_init_seq
diff --git a/red/_validation.py b/red/_validation.py
index 9381ffd..6e5a0e4 100644
--- a/red/_validation.py
+++ b/red/_validation.py
@@ -11,12 +11,11 @@
 def check_data(
     data: _npt.NDArray[_np.float64], one_dim_allowed: bool = False
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Assert that data passed is a numpy array where
-    the first dimension is the number of chains and
-    the second dimension is the number of samples.
-    If the array is one dimensional, add a second
-    dimension with length 1.
+    """Validate and reshape input data to 2D ``(n_chains, n_samples)``.
+
+    Asserts that the data passed is a numpy array where the first dimension is the number
+    of chains and the second dimension is the number of samples. If the array is one
+    dimensional, a second dimension with length 1 is added.
 
     Parameters
     ----------
@@ -30,6 +29,7 @@ def check_data(
     -------
     np.ndarray
         Data with shape (n_chains, n_samples).
+
     """
     # Check that data is a numpy array.
     if not isinstance(data, _np.ndarray):
diff --git a/red/_version.py b/red/_version.py
index 3c07d0f..207b274 100644
--- a/red/_version.py
+++ b/red/_version.py
@@ -1 +1 @@
-__version__ = "0.1.4+6.gf1edc4b.dirty"
+__version__ = "0.1.4+5.g7888c80.dirty"
diff --git a/red/confidence_intervals.py b/red/confidence_intervals.py
index 6a5bbff..d602a32 100644
--- a/red/confidence_intervals.py
+++ b/red/confidence_intervals.py
@@ -16,9 +16,9 @@ def get_conf_int_init_seq(
     min_max_lag_time: int = 3,
     max_max_lag_time: int | None = None,
 ) -> float:
-    """
-    Calculate the confidence interval for the mean of a time
-    series using initial sequence methods. See Geyer, 1992:
+    """Calculate the confidence interval for the mean of a time series.
+
+    Uses initial sequence methods. See Geyer, 1992:
     https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -47,6 +47,7 @@ def get_conf_int_init_seq(
     -------
     float
         The standard error of the mean.
+
     """
     # Get the correlated estimate of the variance.
     var_cor, max_lag, acovf = _get_variance_initial_sequence(
diff --git a/red/equilibration.py b/red/equilibration.py
index 74ec9ee..9d3cc77 100644
--- a/red/equilibration.py
+++ b/red/equilibration.py
@@ -31,13 +31,12 @@ def detect_equilibration_init_seq(
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
     plot_max_lags: bool = True,
 ) -> tuple[float | int, float, float]:
-    r"""
-    Detect the equilibration time of a time series by finding the minimum
-    squared standard error (SSE), or maximum effective sample size (ESS)
-    of the time series, using initial sequence estimators of the variance.
-    This is done by computing the SSE at each time point, discarding all
-    samples before the time point. The index of the time point with
-    the minimum SSE or maximum ESS is taken to be the point of equilibration.
+    r"""Detect the equilibration time of a time series by finding the minimum SSE or maximum ESS.
+
+    The variance is estimated using initial sequence estimators. This is done by computing
+    the squared standard error (SSE) at each time point, discarding all samples before the
+    time point. The index of the time point with the minimum SSE, or maximum effective
+    sample size (ESS), is taken to be the point of equilibration.
 
     Parameters
     ----------
@@ -105,6 +104,7 @@ def detect_equilibration_init_seq(
 
     equil_ess: float
         The effective sample size at the equilibration point.
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -196,13 +196,12 @@ def detect_equilibration_window(
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
     plot_window_size: bool = True,
 ) -> tuple[float | int, float, float]:
-    r"""
-    Detect the equilibration time of a time series by finding the minimum
-    squared standard error (SSE) or maximum effective sample size (ESS)
-    of the time series, using window estimators of the variance. This is
-    done by computing the SSE at each time point, discarding all samples
-    before the time point. The index of the time point with the minimum
-    SSE is taken to be the point of equilibration.
+    r"""Detect the equilibration time of a time series by finding the minimum SSE or maximum ESS.
+
+    The variance is estimated using window estimators. This is done by computing the
+    squared standard error (SSE) at each time point, discarding all samples before the time
+    point. The index of the time point with the minimum SSE, or maximum effective sample
+    size (ESS), is taken to be the point of equilibration.
 
     Parameters
     ----------
@@ -264,6 +263,7 @@ def detect_equilibration_window(
 
     equil_ess: float
         The effective sample size at the equilibration point.
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -343,11 +343,11 @@ def get_paired_t_p_timeseries(
     final_block_size: float = 0.5,
     t_test_sidedness: str = "two-sided",
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.int64 | _np.float64]]:
-    """
-    Get a timeseries of the p-values from a paired t-test on the differences
-    between sample means between intial and final portions of the data. The timeseries
-    is obtained by repeatedly discarding more data from the time series between
-    calculations of the p-value.
+    """Get a timeseries of the p-values from a paired t-test.
+
+    The p-values come from a paired t-test on the differences between sample means between
+    initial and final portions of the data. The timeseries is obtained by repeatedly
+    discarding more data from the time series between calculations of the p-value.
 
     Parameters
     ----------
@@ -385,6 +385,7 @@ def get_paired_t_p_timeseries(
 
     np.ndarray
         The times at which the p-values were calculated.
+
     """
     # Check that the data is valid.
     data = check_data(data, one_dim_allowed=False)
@@ -474,14 +475,14 @@ def detect_equilibration_paired_t_test(
     time_units: str = "ns",
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> _np.int64 | _np.float64:
-    r"""
-    Detect the equilibration time of a time series by performing a paired
-    t-test between initial and final portions of the time series. This is repeated
-    , discarding more data from the time series between repeats. If the p-value
-    is greater than the threshold, there is no significant evidence that the data is
-    no equilibrated and the timeseries is taken to be equilibrated at this time
-    point. This test may be useful when we care only about systematic bias in the
-    data, and do not care about detecting inter-run differences.
+    r"""Detect the equilibration time of a time series by performing a paired t-test.
+
+    A paired t-test is performed between initial and final portions of the time series.
+    This is repeated, discarding more data from the time series between repeats. If the
+    p-value is greater than the threshold, there is no significant evidence that the data
+    is not equilibrated and the timeseries is taken to be equilibrated at this time point.
+    This test may be useful when we care only about systematic bias in the data, and do
+    not care about detecting inter-run differences.
 
     Parameters
     ----------
@@ -533,6 +534,7 @@ def detect_equilibration_paired_t_test(
     np.float64 | np.int64
         The time (or index, if no times are supplied) at which
         the time series is equilibrated.
+
     """
     # Validate data.
     data = check_data(data, one_dim_allowed=False)
diff --git a/red/ess.py b/red/ess.py
index 87beef5..b9dfe94 100644
--- a/red/ess.py
+++ b/red/ess.py
@@ -14,21 +14,21 @@
 def convert_sse_series_to_ess_series(
     data: _npt.NDArray[_np.float64], sse_series: _npt.NDArray[_np.float64]
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Convert a series of squared standard errors to a series of effective sample sizes.
+    """Convert a series of squared standard errors to a series of effective sample sizes.
 
     Parameters
     ----------
+    data : np.ndarray
+        The time series data with shape (n_chains, n_samples).
+
     sse_series : np.ndarray
         The squared standard error series.
 
-    uncor_vars : np.ndarray
-        The uncorrelated variances.
-
     Returns
     -------
     np.ndarray
         The effective sample size series.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
@@ -51,10 +51,10 @@ def get_ess_series_init_seq(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of effective sample sizes for a time series as data
-    is discarded from the beginning of the time series. The autocorrelation
-    is computed using the sequence estimator specified.
+    """Compute a series of effective sample sizes for a time series.
+
+    Effective sample sizes are computed as data is discarded from the beginning of the
+    time series. The autocorrelation is computed using the sequence estimator specified.
 
     Parameters
     ----------
@@ -93,6 +93,7 @@ def get_ess_series_init_seq(
 
     np.ndarray
         The maximum lag times used.
+
     """
     sse_series, max_lag_times = _get_sse_series_init_seq(
         data,
@@ -114,10 +115,11 @@ def get_ess_series_window(
     window_size_fn: _Callable[[int], int] | None = lambda x: round(x**0.5),
     window_size: int | None = None,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of effective sample sizes for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the window size and kernel specified.
+    """Compute a series of effective sample sizes for a time series.
+
+    Effective sample sizes are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the window size and kernel
+    specified.
 
     Parameters
     ----------
@@ -142,6 +144,7 @@ def get_ess_series_window(
 
     np.ndarray
         The window sizes used.
+
     """
     sse_series, max_lag_times = _get_sse_series_window(
         data, kernel=kernel, window_size_fn=window_size_fn, window_size=window_size
@@ -153,10 +156,10 @@ def get_ess_series_window(
 
 
 def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the statistical inefficiency of a time series by dividing
-    the inter-run variance estimate by the intra-run variance estimate.
-    More than one run is required.
+    """Compute the statistical inefficiency of a time series.
+
+    The statistical inefficiency is computed by dividing the inter-run variance estimate
+    by the intra-run variance estimate. More than one run is required.
 
     Parameters
     ----------
@@ -169,6 +172,7 @@ def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) ->
     -------
     float
         The statistical inefficiency.
+
     """
     g = inter_run_variance(data) / intra_run_variance(data)
     # Ensure that the statistical inefficiency is at least 1.
@@ -178,10 +182,11 @@ def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) ->
 def statistical_inefficiency_lugsail_variance(
     data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3
 ) -> float:
-    """
-    Compute the statistical inefficiency of a time series by dividing
-    the lugsail replicated batch means variance estimate by the
-    intra-run variance estimate. This is applicable to a single run.
+    """Compute the statistical inefficiency of a time series.
+
+    The statistical inefficiency is computed by dividing the lugsail replicated batch
+    means variance estimate by the intra-run variance estimate. This is applicable to a
+    single run.
 
     Parameters
     ----------
@@ -198,6 +203,7 @@ def statistical_inefficiency_lugsail_variance(
     -------
     float
         The statistical inefficiency.
+
     """
     g = lugsail_variance(data, n_pow=n_pow) / intra_run_variance(data)
     # Ensure that the statistical inefficiency is at least 1.
@@ -205,11 +211,11 @@ def statistical_inefficiency_lugsail_variance(
 
 
 def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the effective sample size of a time series by dividing
-    the total number of samples by the statistical inefficiency, where
-    the statistical inefficiency is calculated using the ratio of the
-    inter-run and intra-run variance estimates.
+    """Compute the effective sample size of a time series.
+
+    The effective sample size is computed by dividing the total number of samples by the
+    statistical inefficiency, where the statistical inefficiency is calculated using the
+    ratio of the inter-run and intra-run variance estimates.
 
     Parameters
     ----------
@@ -222,6 +228,7 @@ def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The effective sample size.
+
     """
     data = check_data(data, one_dim_allowed=False)
     n_runs: int = data.shape[0]
@@ -231,11 +238,11 @@ def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def ess_lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Compute the effective sample size of a time series by dividing
-    the total number of samples by the statistical inefficiency, where
-    the statistical inefficiency is calculated using the ratio of the
-    lugsail replicated batch means and intra-run variance estimates.
+    """Compute the effective sample size of a time series.
+
+    The effective sample size is computed by dividing the total number of samples by the
+    statistical inefficiency, where the statistical inefficiency is calculated using the
+    ratio of the lugsail replicated batch means and intra-run variance estimates.
 
     Parameters
     ----------
@@ -252,6 +259,7 @@ def ess_lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3)
     -------
     float
         The effective sample size.
+
     """
     data = check_data(data, one_dim_allowed=True)
     n_runs: int = data.shape[0]
diff --git a/red/gelman_rubin.py b/red/gelman_rubin.py
index 03f4480..6bae94b 100644
--- a/red/gelman_rubin.py
+++ b/red/gelman_rubin.py
@@ -16,10 +16,10 @@
 
 
 def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the Gelman-Rubin diagnostic according to
-    equation 4 in  Statist. Sci. 36(4): 518-529
-    (November 2021). DOI: 10.1214/20-STS812
+    """Compute the Gelman-Rubin diagnostic.
+
+    Uses equation 4 in Statist. Sci. 36(4): 518-529 (November 2021).
+    DOI: 10.1214/20-STS812.
 
     Parameters
     ----------
@@ -32,6 +32,7 @@ def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The Gelman-Rubin diagnostic.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=False)
@@ -53,11 +54,10 @@ def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def stable_gelman_rubin(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Compute the stable Gelman-Rubin diagnostic according to
-    equation 7 in  Statist. Sci. 36(4): 518-529
-    (November 2021). DOI: 10.1214/20-STS812. This is applicable to
-    a single run.
+    """Compute the stable Gelman-Rubin diagnostic.
+
+    Uses equation 7 in Statist. Sci. 36(4): 518-529 (November 2021).
+    DOI: 10.1214/20-STS812. This is applicable to a single run.
     """
     # Validate the data.
     data = _check_data(data, one_dim_allowed=True)
diff --git a/red/plot.py b/red/plot.py
index 567c088..5f34097 100644
--- a/red/plot.py
+++ b/red/plot.py
@@ -24,8 +24,7 @@ def plot_timeseries(
     time_units: str = "ns",
     y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> None:
-    r"""
-    Plot the (multi-run) time series data.
+    r"""Plot the (multi-run) time series data.
 
     Parameters
     ----------
@@ -51,6 +50,7 @@ def plot_timeseries(
 
     y_label : str, optional
         The y-axis label. The default is "$\Delta G$ / kcal mol$^{-1}$".
+
     """
     # Check that data is valid.
     data = check_data(data, one_dim_allowed=True)
@@ -118,8 +118,7 @@ def plot_p_values(
     time_units: str = "ns",
     threshold_times: _npt.NDArray[_np.int64 | _np.float64] | None = None,
 ) -> None:
-    """
-    Plot the p-values of the paired t-test.
+    """Plot the p-values of the paired t-test.
 
     Parameters
     ----------
@@ -142,6 +141,7 @@ def plot_p_values(
         The times to plot the p-value threshold at. If None, this is
         set to the times at which the data was sampled. Useful when
         using this plot underneath a time series plot.
+
     """
     # Check that p_values is valid.
     if not isinstance(p_values, _np.ndarray) or not isinstance(times, _np.ndarray):
@@ -204,8 +204,7 @@ def plot_sse(
     variance_y_label: str = r"$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$",
     reciprocal: bool = True,
 ) -> tuple[list[_Artist], list[_Any]]:
-    r"""
-    Plot the squared standard error (SSE) estimate against time.
+    r"""Plot the squared standard error (SSE) estimate against time.
 
     Parameters
     ----------
@@ -229,7 +228,7 @@ def plot_sse(
 
     variance_y_label : str, optional
         The y-axis label for the variance. The default is
-        "$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$".
+        "$\frac{1}{\\sigma^2(\\Delta G)}$ / kcal$^{-2}$ mol$^2$".
 
     reciprocal : bool, optional, default=True
         Whether to plot the reciprocal of the SSE.
@@ -241,6 +240,7 @@ def plot_sse(
 
     labels : List[str]
         The labels for the legend.
+
     """
     # Check that sse is valid.
     if not isinstance(sse, _np.ndarray) or not isinstance(times, _np.ndarray):
@@ -315,16 +315,16 @@ def plot_equilibration_paired_t_test(
     time_units: str = "ns",
     data_y_label: str = r"$\Delta G$ / kcal mol$^{-1}$",
 ) -> tuple[_Axes, _Axes]:
-    r"""
-    Plot the p-values of the paired t-test against time, underneath the
-    time series data.
+    r"""Plot the p-values of the paired t-test against time.
+
+    The p-values are plotted underneath the time series data.
 
     Parameters
     ----------
     fig : plt.Figure
         The figure to plot on.
 
-    gridspec_obj : plt.GridSpec
+    subplot_spec : plt.GridSpec
         The gridspec to use for the plot.
 
     data : np.ndarray
@@ -358,6 +358,7 @@ def plot_equilibration_paired_t_test(
 
     ax_bottom : Axes
         The axes for the p-values.
+
     """
     with _plt.style.context(PLT_STYLE):
         # We need to split the gridspec into two subplots, one for the time series data (above)
@@ -415,9 +416,9 @@ def plot_equilibration_min_sse(
     variance_y_label: str = r"$\frac{1}{\sigma^2(\Delta G)}$ / kcal$^{-2}$ mol$^2$",
     reciprocal: bool = True,
 ) -> tuple[_Axes, _Axes]:
-    r"""
-    Plot the (reciprocal of the) squared standard error (SSE)
-    estimates against time, underneath the time series data.
+    r"""Plot the (reciprocal of the) squared standard error (SSE) estimates against time.
+
+    The estimates are plotted underneath the time series data.
 
     Parameters
     ----------
@@ -472,6 +473,7 @@ def plot_equilibration_min_sse(
 
     ax_bottom : Axes
         The axes for the p-values.
+
     """
     with _plt.style.context(PLT_STYLE):
         data = check_data(data, one_dim_allowed=True)
diff --git a/red/sse.py b/red/sse.py
index 29a9018..ab13033 100644
--- a/red/sse.py
+++ b/red/sse.py
@@ -17,10 +17,11 @@ def get_sse_series_init_seq(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of squared standard errors for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the sequence estimator specified.
+    """Compute a series of squared standard errors for a time series.
+
+    Squared standard errors are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the sequence estimator
+    specified.
 
     Parameters
     ----------
@@ -59,6 +60,7 @@ def get_sse_series_init_seq(
 
     np.ndarray
         The maximum lag times used.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
@@ -89,10 +91,11 @@ def get_sse_series_window(
     window_size: int | None = None,
     frac_padding: float = 0.1,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Compute a series of squared standard errors for a time series as data
-    is discarded from the beginning of the time series. The squared standard
-    error is computed using the window size and kernel specified.
+    """Compute a series of squared standard errors for a time series.
+
+    Squared standard errors are computed as data is discarded from the beginning of the
+    time series. The squared standard error is computed using the window size and kernel
+    specified.
 
     Parameters
     ----------
@@ -123,6 +126,7 @@ def get_sse_series_window(
 
     np.ndarray
         The window sizes used.
+
     """
     # Validate the data.
     data = check_data(data, one_dim_allowed=True)
diff --git a/red/variance.py b/red/variance.py
index 0806b2d..206cbf5 100644
--- a/red/variance.py
+++ b/red/variance.py
@@ -1,5 +1,4 @@
-"""
-Functions to calculate the variance of a time series, accounting for autocorrelation.
+"""Functions to calculate the variance of a time series, accounting for autocorrelation.
 
 Methods implemented:
 
@@ -31,8 +30,8 @@
 def _compute_autocovariance_no_fft(
     data: _npt.NDArray[_np.float64], max_lag: int
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the auto-covariance as a function of lag time for a time series.
+    """Calculate the auto-covariance as a function of lag time for a time series.
+
     Avoids using statsmodel's acovf function as using numpy's dot function and jit
     gives a substantial speedup.
 
@@ -50,6 +49,7 @@ def _compute_autocovariance_no_fft(
     -------
     numpy.ndarray
         The auto-correlation function of the time series.
+
     """
     # Don't use statsmodel's acovf as we can get a substantial speedup by using
     # numpy's dot function and jit.
@@ -70,10 +70,11 @@ def _compute_autocovariance_no_fft(
 def _compute_autocovariance_fft(
     data: _npt.NDArray[_np.float64], max_lag: int
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the autocovariance using the FFT method, as implemented in statsmodels.
+    """Calculate the autocovariance using the FFT method, as implemented in statsmodels.
+
     Note that we can speed this up for large arrays by rewriting to directly use numpy's fft
     function and using jit with rocket-fft https://github.com/styfenschaer/rocket-fft.
+
     Parameters
     ----------
     data : numpy.ndarray
@@ -88,6 +89,7 @@ def _compute_autocovariance_fft(
     -------
     numpy.ndarray
     The auto-correlation function of the time series.
+
     """
     autocov_fn = _acovf(data, adjusted=False, nlag=max_lag, fft=True, demean=False)
     autocov_fn = autocov_fn.astype(_np.float64, copy=False)
@@ -100,8 +102,7 @@ def _get_autocovariance(
     mean: None | float = None,
     fft: bool = False,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the auto-covariance as a function of lag time for a time series.
+    """Calculate the auto-covariance as a function of lag time for a time series.
 
     Parameters
     ----------
@@ -126,6 +127,7 @@ def _get_autocovariance(
     -------
     numpy.ndarray
         The auto-correlation function of the time series.
+
     """
     # Copy the data so we don't modify the original.
     data = data.copy()
@@ -153,8 +155,7 @@ def _get_autocovariance(
 def _get_gamma_cap(
     autocov_series: _npt.NDArray[_np.float64],
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Compute the capitial gamma function from the auto-covariance function.
+    """Compute the capitial gamma function from the auto-covariance function.
 
     Parameters
     ----------
@@ -165,6 +166,7 @@ def _get_gamma_cap(
     -------
     numpy.ndarray
         The capital gamma function of the time series.
+
     """
     # Get the length of the time series.
     n_samples = autocov_series.shape[0]
@@ -191,8 +193,8 @@ def _get_initial_positive_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """ "
-    Get the initial positive sequence from the capital gamma function of a time series.
+    """Get the initial positive sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -208,6 +210,7 @@ def _get_initial_positive_sequence(
     -------
     numpy.ndarray
         The initial positive sequence.
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     # gamma_cap = _deepcopy(gamma_cap)
@@ -227,8 +230,8 @@ def _get_initial_monotone_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Get the initial monotone sequence from the capital gamma function of a time series.
+    """Get the initial monotone sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -243,6 +246,7 @@ def _get_initial_monotone_sequence(
     -------
     numpy.ndarray
         The initial monotone sequence.
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     gamma_cap = gamma_cap.copy()
@@ -263,8 +267,8 @@ def _get_initial_convex_sequence(
     gamma_cap: _npt.NDArray[_np.float64],
     min_max_lag_time: int = 3,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Get the initial convex sequence from the capital gamma function of a time series.
+    """Get the initial convex sequence from the capital gamma function of a time series.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -286,6 +290,7 @@ def _get_initial_convex_sequence(
     MIT License.
     YEAR: 2005, 2009, 2010, 2012
     COPYRIGHT HOLDER: Charles J. Geyer and Leif T. Johnson
+
     """
     # Make a copy of gamma_cap so we don't modify the original.
     gamma_con = gamma_cap.copy()
@@ -342,8 +347,7 @@ def _get_autocovariance_window(
     kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,
     window_size: int = 10,
 ) -> _npt.NDArray[_np.float64]:
-    """
-    Calculate the autocovariance of a time series using window estimators.
+    """Calculate the autocovariance of a time series using window estimators.
 
     Parameters
     ----------
@@ -362,6 +366,7 @@ def _get_autocovariance_window(
     -------
     np.ndarray
         The autocovariance of the time series as a function of lag time.
+
     """
     n_runs, n_samples = data.shape
     if n_samples < window_size:
@@ -391,10 +396,11 @@ def _get_autocovariance_window(
 def _smoothen_max_lag_times(
     max_lag_times: _npt.NDArray[_np.float64],
 ) -> _npt.NDArray[_np.int64]:
-    """
-    Smoothen a list of maximum lag times by a) converting them to a monotinically
-    decreasing sequence and b) linearly interpolating between points where the sequence
-    changes. This may be useful when the max lag times are noisy.
+    """Smoothen a list of maximum lag times.
+
+    This is done by a) converting them to a monotonically decreasing sequence and b)
+    linearly interpolating between points where the sequence changes. This may be useful
+    when the max lag times are noisy.
 
     Parameters
     ----------
@@ -405,6 +411,7 @@ def _smoothen_max_lag_times(
     -------
     numpy.ndarray
         The smoothened maximum lag times.
+
     """
     # Get a monotinically decreasing sequence.
     max_lag_times_monotonic = _get_initial_monotone_sequence(max_lag_times, min_max_lag_time=0)
@@ -439,8 +446,8 @@ def get_variance_initial_sequence(
     max_max_lag_time: int | None = None,
     autocov: _npt.NDArray[_np.float64] | None = None,
 ) -> tuple[float, int, _npt.NDArray[_np.float64]]:
-    """
-    Calculate the variance of a time series using initial sequence methods.
+    """Calculate the variance of a time series using initial sequence methods.
+
     See Geyer, 1992: https://www.jstor.org/stable/2246094.
 
     Parameters
@@ -476,6 +483,7 @@ def get_variance_initial_sequence(
 
     numpy.ndarray
         The auto-covariance function of the time series.
+
     """
     # Validate the data.
     data = _check_data(data, one_dim_allowed=True)
@@ -605,10 +613,11 @@ def get_variance_series_initial_sequence(
     smooth_lag_times: bool = False,
     frac_padding: float = 0.1,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Repeatedly calculate the variance of a time series while discarding increasing
-    numbers of samples from the start of the time series. The variance is calculated
-    using initial sequence methods. See Geyer, 1992: https://www.jstor.org/stable/2246094.
+    """Repeatedly calculate the variance of a time series while discarding samples.
+
+    Increasing numbers of samples are discarded from the start of the time series. The
+    variance is calculated using initial sequence methods. See Geyer, 1992:
+    https://www.jstor.org/stable/2246094.
 
     Parameters
     ----------
@@ -647,6 +656,7 @@ def get_variance_series_initial_sequence(
 
     numpy.ndarray
         The maximum lag time used when calculating the auto-correlated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -734,8 +744,7 @@ def get_variance_window(
     kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,
     window_size: int = 10,
 ) -> float:
-    """
-    Calculate the variance of a time series using window estimators.
+    """Calculate the variance of a time series using window estimators.
 
     Parameters
     ----------
@@ -753,6 +762,7 @@ def get_variance_window(
     -------
     float
         The estimated variance of the time series.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -796,10 +806,10 @@ def get_variance_series_window(
     window_size: int | None = None,
     frac_padding: float = 0.1,
 ) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
-    """
-    Repeatedly calculate the variance of a time series while discarding increasing
-    numbers of samples from the start of the time series. The variance is calculated
-    using window estimators.
+    """Repeatedly calculate the variance of a time series while discarding samples.
+
+    Increasing numbers of samples are discarded from the start of the time series. The
+    variance is calculated using window estimators.
 
     Parameters
     ----------
@@ -830,6 +840,7 @@ def get_variance_series_window(
 
     numpy.ndarray
         The window size used at each index.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -882,8 +893,8 @@ def get_variance_series_window(
 
 
 def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size: int) -> float:
-    """
-    Estimate the variance of a time series using the replicated batch means method.
+    """Estimate the variance of a time series using the replicated batch means method.
+
     See section 3.1 in Statist. Sci. 36(4): 518-529 (November 2021).
     DOI: 10.1214/20-STS812 .
 
@@ -899,6 +910,7 @@ def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size:
     -------
     float
         The estimated variance.
+
     """
     data = _check_data(data, one_dim_allowed=True)
 
@@ -929,8 +941,8 @@ def replicated_batch_means_variance(data: _npt.NDArray[_np.float64], batch_size:
 
 
 def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
-    """
-    Estimate the variance of a time series using the lugsail method.
+    """Estimate the variance of a time series using the lugsail method.
+
     See section 3.2 in Statist. Sci. 36(4): 518-529 (November 2021).
     DOI: 10.1214/20-STS812 .
 
@@ -947,6 +959,7 @@ def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> f
     -------
     float
         The estimated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)
@@ -977,9 +990,7 @@ def lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> f
 
 
 def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the variance based on the inter-run differences
-    between means.
+    """Compute the variance based on the inter-run differences between means.
 
     Parameters
     ----------
@@ -990,6 +1001,7 @@ def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The estimated variance.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=False)
@@ -1005,8 +1017,7 @@ def inter_run_variance(data: _npt.NDArray[_np.float64]) -> float:
 
 
 def intra_run_variance(data: _npt.NDArray[_np.float64]) -> float:
-    """
-    Compute the average intra-run variance estimate.
+    """Compute the average intra-run variance estimate.
 
     Parameters
     ----------
@@ -1017,6 +1028,7 @@ def intra_run_variance(data: _npt.NDArray[_np.float64]) -> float:
     -------
     float
         The mean intra-run variance estimate.
+
     """
     # Check that the data is valid.
     data = _check_data(data, one_dim_allowed=True)