diff --git a/.github/workflows/pytest.yaml b/.github/workflows/pytest.yaml
index 67a648e..34d6444 100644
--- a/.github/workflows/pytest.yaml
+++ b/.github/workflows/pytest.yaml
@@ -43,4 +43,4 @@ jobs:
       - name: Install package with test dependencies
         run: pip install -e .[tests]
       - name: Run pytest
-        run: pytest -vv
+        run: pytest -vv ${{ matrix.python-version == '3.10' && '-W ignore::Warning' || '' }}
diff --git a/docs/clean.sh b/docs/clean.sh
index d53c09f..67c7a23 100755
--- a/docs/clean.sh
+++ b/docs/clean.sh
@@ -1,2 +1,2 @@
 #!/usr/bin/env bash
-rm -rfv build
+rm -rfv build jupyter_execute
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 27a6162..64a9731 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -190,6 +190,7 @@ def run(self, **kwargs):
 add_module_names = False
 autodoc_default_options = {
     "undoc-members": True,
+    "special-members": "__call__",
     "members": None,
     "ignore-module-all": True,
 }
diff --git a/src/stacie/cost.py b/src/stacie/cost.py
index aa33cb9..8db015c 100644
--- a/src/stacie/cost.py
+++ b/src/stacie/cost.py
@@ -47,7 +47,7 @@ class LowFreqCost:
     model: SpectrumModel = attrs.field()
     """The model to be fitted to the spectrum."""
 
-    def __call__(self, pars: NDArray[float], *, deriv: int = 0) -> float:
+    def __call__(self, pars: NDArray[float], *, deriv: int = 0) -> list[NDArray[float]]:
         """Evaluate the cost function and its derivatives.
 
         Parameters
@@ -143,7 +143,7 @@ def expected(self, pars: NDArray[float]) -> NDArray[float]:
 
 def logpdf_gamma(
     x: NDArray[float], alpha: NDArray[float], theta: NDArray[float], *, deriv: int = 0
-):
+) -> list[NDArray[float]]:
     """Compute the logarithm of the probability density function of the Gamma distribution.
 
     Parameters
@@ -179,7 +179,9 @@ def logpdf_gamma(
     return results
 
 
-def entropy_gamma(alpha: NDArray[float], theta: NDArray[float], *, deriv: int = 0):
+def entropy_gamma(
+    alpha: NDArray[float], theta: NDArray[float], *, deriv: int = 0
+) -> list[NDArray[float]]:
     """Compute the entropy of the Gamma distribution.
 
     Parameters
diff --git a/src/stacie/cutoff.py b/src/stacie/cutoff.py
index b9fb212..6137df3 100644
--- a/src/stacie/cutoff.py
+++ b/src/stacie/cutoff.py
@@ -61,7 +61,7 @@ def __call__(
         self,
         spectrum: Spectrum,
         model: SpectrumModel,
-        props: dict[str, NDArray],
+        props: dict[str, NDArray[float]],
     ) -> dict[str, float]:
         """Compute the criterion for the given spectrum and model.
 
@@ -72,7 +72,8 @@ def __call__(
         model
             The model to be fitted to the spectrum.
         props
-            The property dictionary returned by the :py:meth:`stacie.cost.LowFreqCost.props` method.
+            The property dictionary being constructed in the
+            :py:func:`stacie.estimate.fit_model_spectrum` function.
 
         Returns
         -------
@@ -135,7 +136,7 @@ def __call__(
         self,
         spectrum: Spectrum,
         model: SpectrumModel,
-        props: dict[str, NDArray],
+        props: dict[str, NDArray[float]],
     ) -> dict[str, float]:
         """The disparity between fits to two different parts of the spectrum."""
         # Compute weights for the two halves and the model
diff --git a/src/stacie/estimate.py b/src/stacie/estimate.py
index 2adc45f..e758bc3 100644
--- a/src/stacie/estimate.py
+++ b/src/stacie/estimate.py
@@ -49,8 +49,43 @@ class Result:
     props: dict[str] = attrs.field()
     """The properties marginalized over the ensemble of cutoff frequencies.
 
-    Properties of this class derive their results from information in this dictionary.
-    See docstring of :func:`fit_model_spectrum` for details.
+    The following properties documented in :func:`fit_model_spectrum` are estimated as
+    weighted averages over the cutoff frequencies:
+
+    - ``amplitudes_model``: model amplitudes at the included frequencies
+    - ``acint``: autocorrelation integral
+    - ``acint_std``: uncertainty of the autocorrelation integral
+    - ``acint_var``: variance of the autocorrelation integral
+    - ``cost_zscore``: z-score of the cost function
+    - ``criterion_zscore``: z-score of the cutoff criterion
+    - ``fcut``: cutoff frequency
+    - ``pars``: model parameters
+    - ``pars_covar``: covariance matrix of the parameters
+
+    Some properties are not averaged over cutoff frequencies:
+
+    - ``ncut``: number of points included in the fit, i.e. with weight larger than WEIGHT_EPS
+    - ``switch_exponent``: exponent used to construct the cutoff
+    - ``weights``: the weights used to combine the spectrum points in the fit
+
+    When using :class:`stacie.model.LorentzModel`, the following properties are added
+    (derived from the marginalized parameters and their covariance):
+
+    - ``pars_lorentz``: Lorentz parameters (converted from the Padé parameters)
+    - ``pars_lorentz_covar``: covariance matrix of the Lorentz parameters
+    - ``corrtime_exp``: exponential correlation time
+    - ``corrtime_exp_var``: variance of the exponential correlation time
+    - ``corrtime_exp_std``: standard error of the exponential correlation time
+    - ``exp_simulation_time``: recommended simulation time based on the exponential correlation time
+    - ``exp_block_time``: recommended block time based on the exponential correlation time
+
+    When using :class:`stacie.model.ExpPolyModel`, the following additional properties are added
+    (derived from the marginalized parameters and their covariance):
+
+    - ``log_acint``: the logarithm of the autocorrelation integral
+    - ``log_acint_var``: variance of the logarithm of the autocorrelation integral
+    - ``log_acint_std``: standard error of the logarithm of the autocorrelation integral
+
     """
 
     history: list[dict[str]] = attrs.field()
@@ -370,56 +405,58 @@ def fit_model_spectrum(
     Returns
     -------
     props
-        A dictionary containing various intermediate results of the cost function calculation,
-        computed for the optimized parameters.
+        A dictionary containing various intermediate results of the cost function calculation.
         See Notes for details.
 
     Notes
     -----
-    The returned dictionary contains the following items if the fit succeeds:
+    The returned dictionary contains at least the following items, irrespective of
+    whether the fit succeeds or fails:
+
+    - ``fcut``: cutoff frequency used
+    - ``ncut``: number of points included in the fit, i.e. with weight larger than WEIGHT_EPS
+    - ``switch_exponent``: exponent used to construct the cutoff
+    - ``neff``: effective number of points used in the fit (sum of weights)
+    - ``pars_init``: initial guess of the parameters
+    - ``criterion``: value of the cutoff criterion, or infinity if the fit fails.
+    - ``msg``: error message, if the fit fails
 
-    - ``acint``: estimate of the autocorrelation integral
-    - ``acint_var``: variance of the estimate of the autocorrelation integral
-    - ``acint_std``: standard error of the estimate of the autocorrelation integral
+    If the fit succeeds, the following additional statistical estimates are also set:
+
+    - ``acint``: autocorrelation integral
+    - ``acint_var``: variance of the autocorrelation integral
+    - ``acint_std``: standard error of the autocorrelation integral
     - ``cost_value``: cost function value
-    - ``cost_grad``: cost Gradient vector (if ``deriv>=1``)
-    - ``cost_hess``: cost Hessian matrix (if ``deriv==2``)
+    - ``cost_grad``: cost gradient vector (if ``deriv >= 1``)
+    - ``cost_hess``: cost Hessian matrix (if ``deriv == 2``)
     - ``cost_hess_scales``: Hessian rescaling vector, see ``robust_posinv``.
     - ``cost_hess_rescaled_evals``: Rescaled Hessian eigenvalues
-    - ``cost_hess_rescaled_evecs``: Rescaled hessian eigenvectors
+    - ``cost_hess_rescaled_evecs``: Rescaled Hessian eigenvectors
     - ``cost_expected``: expected value of the cost function
     - ``cost_var``: expected variance of the cost function
     - ``cost_zscore``: z-score of the cost function
-    - ``switch_exponent``: exponent used to construct the cutoff
-    - ``criterion``: value of the criterion whose minimizer determines the frequency cutoff
-    - ``criterion_expected``: expected value of the criterion
-    - ``criterion_var``: expected variance of the criterion
-    - ``criterion_zscore``: z-score of the criterion
+    - ``criterion_expected``: expected value of the cutoff criterion
+    - ``criterion_var``: expected variance of the cutoff criterion
+    - ``criterion_zscore``: z-score of the cutoff criterion
     - ``ll``: log likelihood
-    - ``pars_init``: initial guess of the parameters
-    - ``pars``: optimized parameters
-    - ``pars_covar``: covariance matrix of the parameters
+    - ``pars``: model parameters
+    - ``pars_covar``: covariance matrix of the model parameters
 
-    The ``LorentzModel`` has the following additional properties:
+    When using :class:`stacie.model.LorentzModel`, the following estimates are added:
 
-    - ``pars_lorentz``: optimized Lorentz parameters (converted from the Padé parameters)
+    - ``pars_lorentz``: Lorentz parameters (converted from the Padé parameters)
     - ``pars_lorentz_covar``: covariance matrix of the Lorentz parameters
     - ``corrtime_exp``: exponential correlation time, the slowest time scale in the sequences
-    - ``corrtime_exp_var``: variance of the estimated exponential correlation time
-    - ``corrtime_exp_std``: standard error of the estimated exponential correlation time
+    - ``corrtime_exp_var``: variance of the exponential correlation time
+    - ``corrtime_exp_std``: standard error of the exponential correlation time
     - ``exp_simulation_time``: recommended simulation time based on the exponential correlation time
     - ``exp_block_time``: recommended block time based on the exponential correlation time
 
-    The ``ExpPolyModel`` has the following additional properties:
+    When using :class:`stacie.model.ExpPolyModel`, the following estimates are added:
 
     - ``log_acint``: the logarithm of the autocorrelation integral
     - ``log_acint_var``: variance of the logarithm of the autocorrelation integral
     - ``log_acint_std``: standard error of the logarithm of the autocorrelation integral
-
-    If the fit fails, the following properties are set:
-
-    - ``criterion``: infinity
-    - ``msg``: error message
     """
     # Create a switching function for a smooth cutoff
     weights = switch_func(spectrum.freqs, fcut, switch_exponent)
diff --git a/src/stacie/model.py b/src/stacie/model.py
index 09ca37e..9359fd1 100644
--- a/src/stacie/model.py
+++ b/src/stacie/model.py
@@ -31,7 +31,7 @@
 class SpectrumModel:
     """Abstract base class for spectrum models.
 
-    Subclasses must override all methods that raise ``NotImplementedError``.
+    Subclasses must override all methods that raise :class:`NotImplementedError`.
 
     The first parameter must have a property that is used when constructing an initial guess:
     When the first parameter increases, the model should increase everywhere,
diff --git a/src/stacie/spectrum.py b/src/stacie/spectrum.py
index d80c448..cc62177 100644
--- a/src/stacie/spectrum.py
+++ b/src/stacie/spectrum.py
@@ -99,7 +99,7 @@ def compute_spectrum(
     timestep: float = 1,
     include_zero_freq: bool = True,
 ) -> Spectrum:
-    r"""Compute a spectrum and store all inputs for ``estimate_acint`` in a ``Spectrum`` instance.
+    r"""Compute a spectrum and return it as a :class:`Spectrum` object.
 
     The spectrum amplitudes are computed as follows:
 
@@ -165,8 +165,9 @@ def compute_spectrum(
     Returns
     -------
     spectrum
-        A ``Spectrum`` object holding all the inputs needed to estimate
+        A :class:`Spectrum` object holding all the inputs needed to estimate
         the integral of the autocorrelation function.
+        This can be used as input to :func:`stacie.estimate.estimate_acint`.
     """
     # Handle tuple (prefactor, sequences) case
     if isinstance(sequences, tuple):
diff --git a/src/stacie/utils.py b/src/stacie/utils.py
index 58273bf..809f563 100644
--- a/src/stacie/utils.py
+++ b/src/stacie/utils.py
@@ -49,8 +49,8 @@ class UnitConfig:
     - The ``*_unit`` attributes are assumed to have the value of a "display unit" in
       the same internal units.
 
-    For example, if your internal time unit is 1 ps and you want to display times
-    in ns, set ``time_unit=1000.0``, because your display unit (1 ns) is 1000 internal units (1 ps).
+    For example, if your internal time unit is 1 ps and you want times to be reported in ns,
+    set ``time_unit = 1000.0``, because your display unit (1 ns) is 1000 internal units (1 ps).
 
     To make these conventions easy to follow (and to avoid unit hell in general),
     it is recommended to pick consistent internal units for your system.
@@ -61,7 +61,7 @@ class UnitConfig:
       which is also how this class is used in STACIE.
 
     For example, when all variables are in SI base units and you want to display time in ns,
-    frequency in THz, and autocorrelation integrals in cm^2/s, then create a ``UnitConfig``
+    frequency in THz, and autocorrelation integrals in cm^2/s, then create a :class:`UnitConfig`
     as follows:
 
     .. code-block:: python
@@ -205,9 +205,9 @@ class PositiveDefiniteError(ValueError):
 def robust_posinv(matrix: NDArray[float]) -> tuple[NDArray, NDArray, NDArray, NDArray]:
     """Compute the eigenvalues, eigenvectors and inverse of a positive definite symmetric matrix.
 
-    This function is a robust version of ``numpy.linalg.eigh`` and ``numpy.linalg.inv``
+    This function is a robust replacement for :func:`numpy.linalg.eigh` and :func:`numpy.linalg.inv`
     that can handle large variations in order of magnitude of the diagonal elements.
-    If the matrix is not positive definite, a ``ValueError`` is raised.
+    If the matrix is not positive definite, a :class:`ValueError` is raised.
 
     Parameters
     ----------
@@ -249,7 +249,7 @@ def robust_posinv(matrix: NDArray[float]) -> tuple[NDArray, NDArray, NDArray, ND
 def robust_dot(scales, evals, evecs, other):
     """Compute the dot product of a robustly diagonalized matrix with another matrix.
 
-    - The first three arguments are the output of ``robust_posinv``.
+    - The first three arguments are the output of :func:`robust_posinv`.
     - To multiply with the inverse, just use element-wise inversion of ``scales`` and ``evals``.
 
     Parameters