diff --git a/docs/interpolation.rst b/docs/interpolation.rst
index 2a356c538..6bed34075 100644
--- a/docs/interpolation.rst
+++ b/docs/interpolation.rst
@@ -1,22 +1,259 @@
 .. _interpolation:
 
-Interpolation of IRFs
-=====================
+Interpolation and Extrapolation of IRFs
+=======================================
 
-This module contains functions to interpolate from a set of IRFs for different
-conditions to a new IRF. Implementations of interpolation algorithms exist as interpolator
-classes and are applied by top-level estimator classes to IRF components. 
-Direct usage of the interpolator classes is discuraged, as only the estimator classes 
-check the data for consistency.
+.. currentmodule:: pyirf.interpolation 
 
-This can e.g. be used to interpolate IRFs for zenith angles of 20° and 40°
-to 30°.
+This module contains functions to inter- or extrapolate from a set of IRFs for different
+conditions to a new IRF. Implementations of interpolation and extrapolation algorithms 
+exist as interpolator and extrapolator classes and are applied by top-level estimator 
+classes to IRF components. 
+Direct usage of the inter- and extrapolator classes is discouraged, as only the estimator classes 
+check the data for consistency.
 
 Most methods support an arbitrary number of interpolation dimensions although it 
-is strongly advised to limit those for resonable results.
+is strongly advised to limit those for reasonable results.
+The herein provided functionalities can e.g. be used to interpolate the IRF
+for a zenith angle of 30° from available IRFs at 20° and 40°.
+
+
+IRF Component Estimator Classes
+-------------------------------
+
+.. autosummary::
+   :nosignatures:
+
+   EffectiveAreaEstimator       Estimate AEFF tables.
+   RadMaxEstimator              Estimate RadMax tables.
+   EnergyDispersionEstimator    Estimate 2D EDISPs.
+   PSFTableEstimator            Estimate PSF tables.
+
+ 
+
+Inter- and Extrapolation Classes
+--------------------------------
+
+This module provides inter- and extrapolation classes that can be 
+plugged into the estimator classes. 
+Not all of these classes support arbitrary grid-dimensions where the grid 
+in this context is the grid of e.g. observation parameters like zenith angle and 
+magnetic field inclination (this would be a 2D grid) on which template IRFs exist 
+and are meant to be inter- or extrapolated.
+
+For parametrized components (Effective Areas and Rad-Max tables) these classes are:
+
+=============================================    ==================  ============    ==================================================================================================
+**Name**                                         **Type**            **Grid-Dim**    **Note**
+=============================================    ==================  ============    ==================================================================================================
+:any:`GridDataInterpolator`                      Interpolation       Arbitrary       See also :any:`scipy.interpolate.griddata`.
+:any:`ParametrizedNearestSimplexExtrapolator`    Extrapolation       1D or 2D        Linear (1D) or baryzentric (2D) extension outside the grid's convex hull from the nearest simplex.
+:any:`ParametrizedNearestNeighborSearcher`       Nearest Neighbor    Arbitrary       Nearest neighbor finder usable instead of inter- and/or extrapolation. 
+=============================================    ==================  ============    ==================================================================================================
+
+For components represented by discretized PDFs (PSF and EDISP tables) these classes are:
+
+=============================================    ==================  ============    ==============================================================================
+**Name**                                         **Type**            **Grid-Dim**    **Note**
+=============================================    ==================  ============    ==============================================================================
+:any:`QuantileInterpolator`                      Interpolation       Arbitrary       Adaption of [Hol+13]_ and [Rea99]_ to discretized PDFs.
+:any:`MomentMorphInterpolator`                   Interpolation       1D or 2D        Adaption of [Baa+15]_ to discretized PDFs.
+:any:`MomentMorphNearestSimplexExtrapolator`     Extrapolation       1D or 2D        Extension of [Baa+15]_ beyond the grid's convex hull from the nearest simplex.
+:any:`DiscretePDFNearestNeighborSearcher`        Nearest Neighbor    Arbitrary       Nearest neighbor finder usable instead of inter- and/or extrapolation. 
+=============================================    ==================  ============    ==============================================================================
+
+.. [Hol+13] B. E. Hollister and A. T. Pang (2013). Interpolation of Non-Gaussian Probability Distributions for Ensemble Visualization.
+    https://engineering.ucsc.edu/sites/default/files/technical-reports/UCSC-SOE-13-13.pdf
+.. [Rea99] A. L. Read (1999). Linear Interpolation of Histograms.
+    Nucl. Instrum. Methods Phys. Res. A 425, 357-360. https://doi.org/10.1016/S0168-9002(98)01347-3
+.. [Baa+15] M. Baak, S. Gadatsch, R. Harrington and W. Verkerke (2015). Interpolation between 
+    multi-dimensional histograms using a new non-linear moment morphing method
+    Nucl. Instrum. Methods Phys. Res. A 771, 39-48. https://doi.org/10.1016/j.nima.2014.10.033
+
+
+Using Estimator Classes
+-----------------------
+
+Usage of the estimator classes is simple. 
+As an example, consider CTA's Prod5 IRFs [CTA+21]_, they can be downloaded manually or by executing 
+``download_irfs.py`` in ``pyirf's`` root directory, which downloads them to ``.../pyirf/irfs/``.
+The estimator classes can simply be used by first creating an instance of the respective class with all 
+relevant information and then using the object's ``__call__`` interface the obtain results for a specific 
+target point.
+As the energy dispersion represents one of the discretized PDF IRF components, one can use the 
+``MomentMorphInterpolator`` for interpolation and the ``DiscretePDFNearestNeighborSearcher`` 
+for extrapolation.
+
+.. code-block:: python
+
+   import numpy as np
+
+   from gammapy.irf import load_irf_dict_from_file
+   from glob import glob
+   from pyirf.interpolation import (
+      EnergyDispersionEstimator, 
+      MomentMorphInterpolator,
+      DiscretePDFNearestNeighborSearcher
+   )
+
+   # Load IRF data, replace path with actual path
+   PROD5_IRF_PATH = "pyirf/irfs/*.fits.gz"
+
+   irfs = [load_irf_dict_from_file(path) for path in sorted(glob(PROD5_IRF_PATH))]
+
+   edisps = np.array([irf["edisp"].quantity for irf in irfs])
+   bin_edges = irfs[0]["edisp"].axes["migra"].edges
+   # IRFs are for zenith distances of 20, 40 and 60 deg 
+   zen_pnt = np.array([[20], [40], [60]])
+
+   # Create estimator instance
+   edisp_estimator = EnergyDispersionEstimator(
+        grid_points=zen_pnt,
+        migra_bins=bin_edges,
+        energy_dispersion=edisps,
+        interpolator_cls=MomentMorphInterpolator,
+        interpolator_kwargs=None,
+        extrapolator_cls=DiscretePDFNearestNeighborSearcher,
+        extrapolator_kwargs=None,
+    )
+
+    # Estimate energy dispersions
+    interpolated_edisp = edisp_estimator(np.array([[30]]))
+    extrapolated_edisp = edisp_estimator(np.array([[10]]))
+
+
+.. [CTA+21] Cherenkov Telescope Array Observatory & Cherenkov Telescope Array Consortium. (2021).
+   CTAO Instrument Response Functions - prod5 version v0.1 (v0.1) [Data set]. Zenodo. 
+   https://doi.org/10.5281/zenodo.5499840
+
+
+Creating new Estimator Classes
+------------------------------
+
+To create a estimator class for an IRF component not yet implemented, one can simply 
+inherit from respective base class.
+There are two, tailored to either parametrized or discrete PDF components.
+
+.. autosummary::
+   :nosignatures:
+
+   ParametrizedComponentEstimator   Parametrized components 
+   DiscretePDFComponentEstimator    Discrete PDF components
+
+Consider an example, where one is interested in an estimator for simple Gaussians.
+As this is already the scope of the ``DiscretePDFComponentEstimator`` base class and 
+for the sake of this demonstration, let the Gaussians come with some 
+units attached that need handling:
+
+.. code-block:: python
+
+   import astropy.units as u
+   from pyirf.interpolation import (DiscretePDFComponentEstimator,
+                                    MomentMorphInterpolator)
+
+   class GaussianEstimatior(DiscretePDFComponentEstimator):
+      @u.quantity_input(gaussians=u.m)
+      def __init__(
+         self,
+         grid_points,
+         bin_edges,
+         gaussians,
+         interpolator_cls=MomentMorphInterpolator,
+         interpolator_kwargs=None,
+         extrapolator_cls=None,
+         extrapolator_kwargs=None,
+      ):
+         if interpolator_kwargs is None:
+            interpolator_kwargs = {}
+
+         if extrapolator_kwargs is None:
+            extrapolator_kwargs = {}
+
+         self.unit = gaussians.unit
+
+         super().__init__(
+            grid_points=grid_points,
+            bin_edges=bin_edges,
+            binned_pdf=gaussians.to_value(u.m),
+            interpolator_cls=interpolator_cls,
+            interpolator_kwargs=interpolator_kwargs,
+            extrapolator_cls=extrapolator_cls,
+            extrapolator_kwargs=extrapolator_kwargs,
+         )
+
+      def __call__(self, target_point):
+         res = super().__call__(target_point)
+
+         # Return result with correct unit
+         return u.Quantity(res, u.m, copy=False).to(self.unit)
+
+This new estimator class can now be used just like any other estimator class already 
+implemented in ``pyirf.interpolation``. 
+While the ``extrapolator_cls`` argument can be empty when creating an instance of 
+``GaussianEstimator``, effectively disabling extrapolation and raising an error in 
+case it would be needed regardless, assume the desired extrapolation method to be 
+``MomentMorphNearestSimplexExtrapolator``:
+
+.. code-block:: python
+
+   import numpy as np
+   from pyirf.interpolation import MomentMorphNearestSimplexExtrapolator
+   from scipy.stats import norm 
+
+   bins = np.linspace(-10, 10, 51)
+   grid = np.array([[1], [2], [3]])
+
+   gaussians = np.array([np.diff(norm(loc=x, scale=1/x).cdf(bins))/np.diff(bins) for x in grid])
+
+   estimator = GaussianEstimatior(
+      grid_points = grid,
+      bin_edges = bins,
+      gaussians = gaussians * u.m,
+      interpolator_cls = MomentMorphInterpolator,
+      extrapolator_cls = MomentMorphNearestSimplexExtrapolator
+   )
+
+This estimator object can now easily be used to estimate Gaussians at arbitrary target points:
+
+.. code-block:: python
+
+   targets = np.array([[0.9], [1.5]])
+
+   results = u.Quantity([estimator(target).squeeze() for target in targets])
+
+
+Helper Classes
+--------------
+
+.. autosummary::
+   :nosignatures:
+
+   PDFNormalization
+
+
+Base Classes
+------------
+
+.. autosummary::
+   :nosignatures:
+
+   BaseComponentEstimator
+   ParametrizedComponentEstimator
+   DiscretePDFComponentEstimator
+   BaseInterpolator
+   ParametrizedInterpolator
+   DiscretePDFInterpolator
+   BaseExtrapolator
+   ParametrizedExtrapolator
+   DiscretePDFExtrapolator
+   BaseNearestNeighborSearcher
 
 
-Reference/API
--------------
+Full API
+--------
 
 .. automodapi:: pyirf.interpolation
+   :no-heading:
+   :no-main-docstr:
+   :inherited-members:
+   :no-inheritance-diagram:
diff --git a/pyirf/interpolation/base_extrapolators.py b/pyirf/interpolation/base_extrapolators.py
index 28b3b38a8..9edd9933b 100644
--- a/pyirf/interpolation/base_extrapolators.py
+++ b/pyirf/interpolation/base_extrapolators.py
@@ -26,7 +26,7 @@ def __init__(self, grid_points):
         self.grid_points = grid_points
         if self.grid_points.ndim == 1:
             self.grid_points = self.grid_points.reshape(*self.grid_points.shape, 1)
-        self.n_points = self.grid_points.shape[0]
+        self.N = self.grid_points.shape[0]
         self.grid_dim = self.grid_points.shape[1]
 
     @abstractmethod
@@ -60,7 +60,7 @@ def __init__(self, grid_points, params):
 
         Parameters
         ----------
-        grid_points, np.ndarray, shape=(n_points, n_dims)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist
         params: np.ndarray, shape=(n_points, ..., n_params)
             Corresponding parameter values at each point in grid_points.
@@ -84,21 +84,23 @@ class DiscretePDFExtrapolator(BaseExtrapolator):
     Derived from pyirf.interpolation.BaseExtrapolator
     """
 
-    def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA):
+    def __init__(
+        self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA
+    ):
         """DiscretePDFExtrapolator
 
         Parameters
         ----------
-        grid_points : np.ndarray, shape=(n_points, n_dims)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist
-        bin_edges : np.ndarray, shape=(n_bins+1)
+        bin_edges: np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        binned_pdf : np.ndarray, shape=(n_points, ..., n_bins)
+        binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, last dimension has to correspond to number of bins for
             the quantity that should be extrapolated (e.g. the Migra axis for EDisp)
-        normalization : PDFNormalization
+        normalization: PDFNormalization
             How the PDF is normalized
 
         Note
diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py
index dc859c868..1d102d933 100644
--- a/pyirf/interpolation/component_estimators.py
+++ b/pyirf/interpolation/component_estimators.py
@@ -9,8 +9,8 @@
 from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator
 from .base_interpolators import (
     DiscretePDFInterpolator,
-    PDFNormalization,
     ParametrizedInterpolator,
+    PDFNormalization,
 )
 from .griddata_interpolator import GridDataInterpolator
 from .quantile_interpolator import QuantileInterpolator
@@ -37,7 +37,7 @@ class BaseComponentEstimator:
     def __init__(self, grid_points):
         """
         Base __init__, doing sanity checks on the grid, building a
-        triangulated version of the grid and intantiating inter- and extrapolator.
+        triangulated version of the grid and instantiating inter- and extrapolator.
 
         Parameters
         ----------
@@ -46,7 +46,7 @@ def __init__(self, grid_points):
         Raises
         ------
         TypeError:
-            When grid_points is not a np.ndarray
+            When grid_points is not a np.ndarray of float compatible values
         TypeError:
             When grid_point has dtype object
         ValueError:
@@ -370,6 +370,8 @@ def __init__(
 
 
 class EffectiveAreaEstimator(ParametrizedComponentEstimator):
+    """Estimator class for effective area tables (AEFF_2D)."""
+
     @u.quantity_input(effective_area=u.m**2, min_effective_area=u.m**2)
     def __init__(
         self,
@@ -382,8 +384,6 @@ def __init__(
         min_effective_area=1 * u.m**2,
     ):
         """
-        Estimator class for effective areas.
-
         Takes a grid of effective areas for a bunch of different parameters
         and inter-/extrapolates (log) effective areas to given value of
         those parameters.
@@ -471,6 +471,8 @@ def __call__(self, target_point):
 
 
 class RadMaxEstimator(ParametrizedComponentEstimator):
+    """Estimator class for rad-max tables (RAD_MAX, RAD_MAX_2D)."""
+
     def __init__(
         self,
         grid_points,
@@ -481,8 +483,6 @@ def __init__(
         extrapolator_kwargs=None,
     ):
         """
-        Estimator class for RAD_MAX tables.
-
         Takes a grid of rad max values for a bunch of different parameters
         and inter-/extrapolates rad max values to given value of those parameters.
 
@@ -545,6 +545,8 @@ def __call__(self, target_point):
 
 
 class EnergyDispersionEstimator(DiscretePDFComponentEstimator):
+    """Estimator class for energy dispersions (EDISP_2D)."""
+
     def __init__(
         self,
         grid_points,
@@ -557,8 +559,6 @@ def __init__(
         axis=-2,
     ):
         """
-        Estimator class for energy dispersions.
-
         Takes a grid of energy dispersions for a bunch of different parameters and
         inter-/extrapolates energy dispersions to given value of those parameters.
 
@@ -631,6 +631,8 @@ def __call__(self, target_point):
 
 
 class PSFTableEstimator(DiscretePDFComponentEstimator):
+    """Estimator class for point spread function tables (PSF_TABLE)."""
+
     @u.quantity_input(psf=u.sr**-1, source_offset_bins=u.deg)
     def __init__(
         self,
@@ -644,8 +646,6 @@ def __init__(
         axis=-1,
     ):
         """
-        Estimator class for point spread functions.
-
         Takes a grid of psfs or a bunch of different parameters and
         inter-/extrapolates psfs to given value of those parameters.
 
@@ -712,7 +712,7 @@ def __init__(
 
     def __call__(self, target_point):
         """
-        Estimating energy dispersions at target_point, inter-/extrapolates as needed and
+        Estimating psf tables at target_point, inter-/extrapolates as needed and
         specified in __init__.
 
         Parameters
@@ -731,4 +731,6 @@ def __call__(self, target_point):
         interpolated_psf_normed = super().__call__(target_point)
 
         # Undo normalisation to get a proper PSF and return
-        return u.Quantity(np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False)
+        return u.Quantity(
+            np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False
+        )
diff --git a/pyirf/interpolation/griddata_interpolator.py b/pyirf/interpolation/griddata_interpolator.py
index 9130a5701..7d50359a9 100644
--- a/pyirf/interpolation/griddata_interpolator.py
+++ b/pyirf/interpolation/griddata_interpolator.py
@@ -9,14 +9,16 @@
 
 
 class GridDataInterpolator(ParametrizedInterpolator):
+    """ "Wrapper arounf scipy.interpolate.griddata."""
+
     def __init__(self, grid_points, params, **griddata_kwargs):
         """Parametrized Interpolator using scipy.interpolate.griddata
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(N, O)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which interpolation templates exist
-        params: np.ndarray, shape=(N, ...)
+        params: np.ndarray, shape=(n_points, ..., n_params)
             Structured array of corresponding parameter values at each
             point in grid_points.
             First dimesion has to correspond to number of grid_points
@@ -48,12 +50,13 @@ def interpolate(self, target_point):
 
         Parameters
         ----------
-        target_point: np.ndarray, shape=(1, O)
+        target_point: np.ndarray, shape=(1, n_dims)
             Target point for interpolation
 
         Returns
         -------
-        interpolant: np.ndarray
+        interpolant: np.ndarray, shape=(1, ..., n_params)
+            Interpolated parameter values
 
         References
         ----------
diff --git a/pyirf/interpolation/moment_morph_interpolator.py b/pyirf/interpolation/moment_morph_interpolator.py
index f561e050d..77c5000b9 100644
--- a/pyirf/interpolation/moment_morph_interpolator.py
+++ b/pyirf/interpolation/moment_morph_interpolator.py
@@ -16,18 +16,18 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization):
 
     Parameters
     ----------
-    bin_edges: np.ndarray, shape=(M+1)
+    bin_edges: np.ndarray, shape=(n_bins+1)
         Array of common bin-edges for binned_pdf
-    binned_pdf: np.ndarray, shape=(N, ..., M)
+    binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
         PDF values from which to compute mean and std
     normalization : PDFNormalization
         How the PDF is normalized
 
     Returns
     -------
-    mean: np.ndarray, shape=(N, ...)
+    mean: np.ndarray, shape=(n_points, ...)
         Estimated mean for each input template
-    std: np.ndarray, shape=(N, ...)
+    std: np.ndarray, shape=(n_points, ...)
         Estimated standard deviation for each input template. Set to width/2 if only one bin in
         the input template is =/= 0
     """
@@ -51,7 +51,7 @@ def _estimate_mean_std(bin_edges, binned_pdf, normalization):
     if np.any(mask):
         width = np.diff(bin_edges)
         # std of a uniform distribution inside the bin
-        uniform_std = np.broadcast_to(np.sqrt(1/12) * width, binned_pdf[mask].shape)
+        uniform_std = np.broadcast_to(np.sqrt(1 / 12) * width, binned_pdf[mask].shape)
         std[mask] = uniform_std[binned_pdf[mask, :] != 0]
 
     return mean, std
@@ -63,17 +63,17 @@ def _lookup(bin_edges, binned_pdf, x):
 
     Parameters
     ----------
-    bin_edges: np.ndarray, shape=(M+1)
+    bin_edges: np.ndarray, shape=(n_bins+1)
         Array of common bin-edges for binned_pdf
-    binned_pdf: np.ndarray, shape=(N, ..., M)
+    binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
         Array of bin-entries, actual
-    x: numpy.ndarray, shape=(N, ..., M)
-        Array of M points for each input template, where the histogram-value (bin-height) should be found
+    x: numpy.ndarray, shape=(n_points, ..., n_bins)
+        Array of n_bins points for each input template, where the histogram-value (bin-height) should be found
 
     Returns
     -------
-    y: numpy.ndarray, shape=(N, ..., M)
-        Array of the bin-heights at the M points x, set to 0 at each point outside the histogram
+    y: numpy.ndarray, shape=(n_points, ..., n_bins)
+        Array of the bin-heights at the n_bins points x, set to 0 at each point outside the histogram
 
     """
     # Find the bin where each point x is located in
@@ -186,18 +186,18 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization):
 
     Parameters
     ----------
-    bin_edges: np.ndarray, shape=(M+1)
+    bin_edges: np.ndarray, shape=(n_bins+1)
         Array of common bin-edges for binned_pdf
-    binned_pdf: np.ndarray, shape=(N, ..., M)
+    binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
         Array of bin-entries, actual
-    coefficients: np.ndarray, shape=(N)
+    coefficients: np.ndarray, shape=(n_points)
         Estimation coefficients for each entry in binned_pdf
     normalization : PDFNormalization
         How the PDF is normalized
 
     Returns
     -------
-    f_new: numpy.ndarray, shape=(1, M)
+    f_new: numpy.ndarray, shape=(1, n_bins)
         Interpolated histogram
 
     References
@@ -216,7 +216,9 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization):
 
     # Estimate mean and std for each input template histogram. First adaption needed to extend
     # the moment morph procedure to histograms
-    mus, sigs = _estimate_mean_std(bin_edges=bin_edges, binned_pdf=binned_pdf, normalization=normalization)
+    mus, sigs = _estimate_mean_std(
+        bin_edges=bin_edges, binned_pdf=binned_pdf, normalization=normalization
+    )
     coefficients = coefficients.reshape(
         binned_pdf.shape[0], *np.ones(mus.ndim - 1, "int")
     )
@@ -259,6 +261,8 @@ def moment_morph_estimation(bin_edges, binned_pdf, coefficients, normalization):
 
 
 class MomentMorphInterpolator(DiscretePDFInterpolator):
+    """Interpolator class providing Moment Morphing to interpolate discretized PDFs."""
+
     def __init__(
         self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA
     ):
@@ -267,11 +271,11 @@ def __init__(
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(N, ...)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which interpolation templates exist. May be one ot two dimensional.
-        bin_edges: np.ndarray, shape=(M+1)
+        bin_edges: np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        binned_pdf: np.ndarray, shape=(N, ..., M)
+        binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points. Interpolation dimension, meaning the
@@ -342,12 +346,12 @@ def interpolate(self, target_point):
 
         Parameters
         ----------
-        target_point: numpy.ndarray
+        target_point: numpy.ndarray, shape=(1, n_dims)
             Value for which the interpolation is performed (target point)
 
         Returns
         -------
-        f_new: numpy.ndarray, shape=(1,...,M,...)
+        f_new: numpy.ndarray, shape=(1,...,n_bins)
             Interpolated and binned pdf
 
         References
diff --git a/pyirf/interpolation/nearest_neighbor_searcher.py b/pyirf/interpolation/nearest_neighbor_searcher.py
index 6ddcf5df7..4543d4c07 100644
--- a/pyirf/interpolation/nearest_neighbor_searcher.py
+++ b/pyirf/interpolation/nearest_neighbor_searcher.py
@@ -71,12 +71,12 @@ def interpolate(self, target_point):
 
         Parameters
         ----------
-        target_point: numpy.ndarray
+        target_point: numpy.ndarray, shape=(1, n_dims)
             Value for which the nearest neighbor should be found (target point)
 
         Returns
         -------
-        content_new: numpy.ndarray, shape=(1,...,M,...)
+        content_new: numpy.ndarray, shape=(1, ...)
             values at nearest neighbor
 
         Note
@@ -100,7 +100,7 @@ def interpolate(self, target_point):
 class DiscretePDFNearestNeighborSearcher(BaseNearestNeighborSearcher):
     """
     Dummy NearestNeighbor approach usable instead of
-    actual Interpolation/Extrapolation.
+    actual interpolation/extrapolation.
     Compatible with discretized PDF IRF component API.
     """
 
@@ -139,7 +139,7 @@ def __init__(self, grid_points, bin_edges, binned_pdf, norm_ord=2):
 class ParametrizedNearestNeighborSearcher(BaseNearestNeighborSearcher):
     """
     Dummy NearestNeighbor approach usable instead of
-    actual Interpolation/Extrapolation
+    actual interpolation/extrapolation
     Compatible with parametrized IRF component API.
     """
 
diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py
index 3dc8c7aa8..a04481be0 100644
--- a/pyirf/interpolation/nearest_simplex_extrapolator.py
+++ b/pyirf/interpolation/nearest_simplex_extrapolator.py
@@ -7,8 +7,11 @@
 import numpy as np
 from scipy.spatial import Delaunay
 
-
-from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator, PDFNormalization
+from .base_extrapolators import (
+    DiscretePDFExtrapolator,
+    ParametrizedExtrapolator,
+    PDFNormalization,
+)
 from .moment_morph_interpolator import (
     barycentric_2D_interpolation_coefficients,
     linesegment_1D_interpolation_coefficients,
@@ -23,17 +26,19 @@
 
 
 class ParametrizedNearestSimplexExtrapolator(ParametrizedExtrapolator):
+    """Extrapolator class extending linear or baryzentric interpolation outside a grid's convex hull."""
+
     def __init__(self, grid_points, params):
         """
-        Extrapolator class using linear extrapolation in one ore two
+        Extrapolator class using linear or baryzentric extrapolation in one ore two
         grid-dimensions.
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(N, ...)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist. May be one ot two dimensional.
             Have to be sorted in accending order for 1D.
-        params: np.ndarray, shape=(N, ...)
+        params: np.ndarray, shape=(n_points, ...)
             Array of corresponding parameter values at each point in grid_points.
             First dimesion has to correspond to number of grid_points
 
@@ -95,7 +100,7 @@ def extrapolate(self, target_point):
 
         Parameters
         ----------
-        target_point: numpy.ndarray
+        target_point: numpy.ndarray, shape=(1, n_dims)
             Value for which the extrapolation is performed (target point)
 
         Returns
@@ -123,7 +128,11 @@ def extrapolate(self, target_point):
 
 
 class MomentMorphNearestSimplexExtrapolator(DiscretePDFExtrapolator):
-    def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA):
+    """Extrapolator class extending moment morphing interpolation outside a grid's convex hull."""
+
+    def __init__(
+        self, grid_points, bin_edges, binned_pdf, normalization=PDFNormalization.AREA
+    ):
         """
         Extrapolator class extending/reusing parts of Moment Morphing
         by allowing for negative extrapolation coefficients computed
@@ -132,12 +141,12 @@ def __init__(self, grid_points, bin_edges, binned_pdf, normalization=PDFNormaliz
 
         Parameters
         ----------
-        grid_points: np.ndarray, shape=(N, ...)
+        grid_points: np.ndarray, shape=(n_points, n_dims)
             Grid points at which templates exist. May be one ot two dimensional.
             Have to be sorted in accending order for 1D.
-        bin_edges: np.ndarray, shape=(M+1)
+        bin_edges: np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        binned_pdf: np.ndarray, shape=(N, ..., M)
+        binned_pdf: np.ndarray, shape=(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points. Extrapolation dimension, meaning the
@@ -201,12 +210,12 @@ def extrapolate(self, target_point):
 
         Parameters
         ----------
-        target_point: numpy.ndarray
+        target_point: numpy.ndarray, shape=(1, n_dims)
             Value for which the extrapolation is performed (target point)
 
         Returns
         -------
-        values: numpy.ndarray, shape=(1, ..., M)
+        values: numpy.ndarray, shape=(1, ..., n_bins)
             Extrapolated discretized PDFs
 
         """
diff --git a/pyirf/interpolation/quantile_interpolator.py b/pyirf/interpolation/quantile_interpolator.py
index 890624191..deaeae184 100644
--- a/pyirf/interpolation/quantile_interpolator.py
+++ b/pyirf/interpolation/quantile_interpolator.py
@@ -27,18 +27,18 @@ def ppf_values(bin_mids, cdfs, quantiles):
 
     Parameters
     ----------
-    bin_mids: numpy.ndarray, shape=(M)
+    bin_mids: numpy.ndarray, shape=(n_bins)
         Bin-mids for each bin along interpolation axis
 
-    cdfs: numpy.ndarray, shape=(N,...,M)
+    cdfs: numpy.ndarray, shape=(n_points,...,n_bins)
         Corresponding cdf-values for all quantiles
 
-    quantiles: numpy.ndarray, shape=(L)
+    quantiles: numpy.ndarray, shape=(n_quantiles)
         Quantiles for which the ppf-values should be estimated
 
     Returns
     -------
-    ppfs: numpy.ndarray, shape=(1,...,L)
+    ppfs: numpy.ndarray, shape=(1,...,n_quantiles)
         Corresponding ppf-values for all quantiles at the target interpolation point
     """
 
@@ -93,20 +93,20 @@ def pdf_from_ppf(bin_edges, interp_ppfs, quantiles):
 
     Parameters
     ----------
-    bin_edges: numpy.ndarray, shape=(M+1)
+    bin_edges: numpy.ndarray, shape=(n_bins+1)
         Edges of the bins in which the final pdf should be binned
 
-    interp_ppfs: numpy.ndarray, shape=(1,...,L)
+    interp_ppfs: numpy.ndarray, shape=(1,...,n_quantiles)
         Corresponding ppf-values for all quantiles at the target_point,
         not to be confused with QunatileInterpolators self.ppfs, the ppfs
         computed from the input distributions.
 
-    quantiles: numpy.ndarray, shape=(L)
+    quantiles: numpy.ndarray, shape=(n_quantiles)
         Quantiles corresponding to the ppf-values in interp_ppfs
 
     Returns
     -------
-    pdf_values: numpy.ndarray, shape=(1,...,M)
+    pdf_values: numpy.ndarray, shape=(1,...,n_bins)
         Recomputed, binned pdf at target_point
     """
     # recalculate pdf values through numerical differentiation
@@ -152,6 +152,8 @@ def norm_pdf(pdf_values, bin_edges, normalization):
 
 
 class QuantileInterpolator(DiscretePDFInterpolator):
+    """Interpolator class providing quantile interpoalation."""
+
     def __init__(
         self,
         grid_points,
@@ -160,15 +162,14 @@ def __init__(
         quantile_resolution=1e-3,
         normalization=PDFNormalization.AREA,
     ):
-        """BinnedInterpolator constructor
-
+        """
         Parameters
         ----------
-        grid_points : np.ndarray
+        grid_points : np.ndarray, shape=(n_points, n_dims)
             Grid points at which interpolation templates exist
-        bin_edges : np.ndarray
+        bin_edges : np.ndarray, shape=(n_bins+1)
             Edges of the data binning
-        binned_pdf : np.ndarray
+        binned_pdf : np.ndarray, shape(n_points, ..., n_bins)
             Content of each bin in bin_edges for
             each point in grid_points. First dimesion has to correspond to number
             of grid_points, the last axis has to correspond to number
@@ -232,12 +233,12 @@ def interpolate(self, target_point):
 
         Parameters
         ----------
-        target_point: numpy.ndarray, shape=(O)
+        target_point: numpy.ndarray, shape=(1, n_dims)
             Value for which the interpolation is performed (target point)
 
         Returns
         -------
-        f_new: numpy.ndarray, shape=(1,...,M,...)
+        f_new: numpy.ndarray, shape=(1, ..., n_bins)
             Interpolated and binned pdf
 
         References
@@ -257,7 +258,9 @@ def interpolate(self, target_point):
         )[..., 1:]
 
         # Renormalize pdf to sum of 1
-        normed_interpolated_pdfs = norm_pdf(interpolated_pdfs, self.bin_edges[1:], self.normalization)
+        normed_interpolated_pdfs = norm_pdf(
+            interpolated_pdfs, self.bin_edges[1:], self.normalization
+        )
 
         # Re-swap axes and set all nans to zero
         return np.nan_to_num(normed_interpolated_pdfs).reshape(1, *self.input_shape[1:])
diff --git a/pyirf/interpolation/utils.py b/pyirf/interpolation/utils.py
index 5936ec63a..5792a60c8 100644
--- a/pyirf/interpolation/utils.py
+++ b/pyirf/interpolation/utils.py
@@ -22,9 +22,9 @@ def plumb_point_dist(line, target):
 
     Parameters
     ----------
-    line: np.ndarray, shape=(2, M)
-        Array of two points spanning a line segment. Might be in two or three dims M.
-    target: np.ndarray, shape=(M)
+    line: np.ndarray, shape=(2, n_dims)
+        Array of two points spanning a line segment. Might be in two or three dims n_dims.
+    target: np.ndarray, shape=(n_dims)
         Target point, of which the minimal distance to line segement is needed
 
     Returns
@@ -75,9 +75,9 @@ def point_facet_angle(line, target):
 
     Parameters
     ----------
-    line: np.ndarray, shape=(2, M)
-        Array of two points spanning a line segment. Might be in two or three dims M.
-    target: np.ndarray, shape=(M)
+    line: np.ndarray, shape=(2, n_dims)
+        Array of two points spanning a line segment. Might be in two or three dims n_dims.
+    target: np.ndarray, shape=(n_dims)
         Target point, of which the angle is needed.
 
     Returns