Repository: scipy/scipy
Pull request: MAINT: stats: shift `NumericalInverseHermite` to UNU.RAN
#### Reference issue

#14600 
Closes gh-14875

#### What does this implement/fix?

Begin translation of FNI to UNU.RAN. Minimal changes have been done to the existing docs, tests, and benchmarks. What's left:

- [x] Tutorial
- [x] Stronger tests
- [x] More docs
- [x] Doc deprecate `tol` parameter
- [x] Add type hints

Two tests with SciPy's distribution suite are failing: one with round-off numerical errors warning (`kappa4` distribution) and the other with an error that CDF was less than 0 (`vonmises` distribution).

A `get_qrng` method has been added to the `_URNG` class to add support for the `qrvs` method of `NumericalInverseHermite`. It samples an array of quasi-uniform random numbers and stores them in a  thread-safe array which is accessed by UNU.RAN while transforming those samples to actual samples using the approximated PPF. As UNU.RAN doesn't provide the truncated boundaries of the distribution used by the generator, we can use the `ppf` method directly (as previously done) to transform the samples.

Backward incompatible changes:

- An array of midpoint u-error was previously calculated and stored in the attribute `fni.midpoint_error`. But there is no way to access the midpoint u-error in UNU.RAN so the attribute has been removed. A `u_error` method has been added instead. Cross-ref gh-14875

#### Additional information
<!--Any additional information you think is important.-->

Problem statement:
BUG: midpoint_error of NumericalInverseHermite returns array instead of float
### Describe your issue.

The midpoint_error attribute should return the maximum error (float), not the array of all errors.

### Reproducing Code Example

```python
from scipy import stats

dist = stats.norm()
rng = stats.NumericalInverseHermite(dist)
rng.midpoint_error
```


### Error message

```shell
No error msg...
```


### SciPy/NumPy/Python version information

SciPy 1.7.0

Hints from issue discussion:

Commit message: DOC: stats: add a note that tol has been doc-deprecated

Patch to review:
diff --git a/benchmarks/benchmarks/stats.py b/benchmarks/benchmarks/stats.py
index aa185a13cfa5..0f105b3948ee 100644
--- a/benchmarks/benchmarks/stats.py
+++ b/benchmarks/benchmarks/stats.py
@@ -494,32 +494,6 @@ def time_halton(self, d, scramble, n, workers):
         seq.random(n, workers=workers)
 
 
-class NumericalInverseHermite(Benchmark):
-
-    param_names = ['distribution']
-    params = [distcont]
-
-    def setup(self, *args):
-        self.rand = [np.random.normal(loc=i, size=1000) for i in range(3)]
-
-    def time_fni(self, distcase):
-        distname, shapes = distcase
-        slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
-        fail_dists = {'beta', 'gausshyper', 'geninvgauss', 'ncf', 'nct',
-                      'norminvgauss', 'genhyperbolic', 'studentized_range'}
-
-        if distname in slow_dists or distname in fail_dists:
-            raise NotImplementedError("skipped")
-
-        dist = getattr(stats, distname)(*shapes)
-
-        with np.testing.suppress_warnings() as sup:
-            sup.filter(RuntimeWarning, "overflow encountered")
-            sup.filter(RuntimeWarning, "divide by zero")
-            sup.filter(RuntimeWarning, "invalid value encountered")
-            stats.NumericalInverseHermite(dist)
-
-
 class DistanceFunctions(Benchmark):
     param_names = ['n_size']
     params = [
diff --git a/benchmarks/benchmarks/stats_sampling.py b/benchmarks/benchmarks/stats_sampling.py
index 4367a8b44d85..d22efa5ecf23 100644
--- a/benchmarks/benchmarks/stats_sampling.py
+++ b/benchmarks/benchmarks/stats_sampling.py
@@ -10,21 +10,16 @@
 # Beta distribution with a = 2, b = 3
 class contdist1:
     def pdf(self, x):
-        if 0 < x < 1:
-            return x * (1-x)**2
-        return 0
+        return 12 * x * (1-x)**2
 
     def dpdf(self, x):
-        if 0 < x < 1:
-            return (1-x)**2 - 2*x*(1-x)
-        return 0
+        return 12 * ((1-x)**2 - 2*x*(1-x))
 
     def cdf(self, x):
-        if x < 0:
-            return 0
-        if x > 1:
-            return 1
-        return stats.beta._cdf(x, 2, 3)
+        return 12 * (x**2/2 - x**3/3 + x**4/4)
+
+    def support(self):
+        return 0, 1
 
     def __repr__(self):
         # asv prints this.
@@ -34,10 +29,10 @@ def __repr__(self):
 # Standard Normal Distribution
 class contdist2:
     def pdf(self, x):
-        return np.exp(-0.5 * x*x)
+        return 1./np.sqrt(2*np.pi) * np.exp(-0.5 * x*x)
 
     def dpdf(self, x):
-        return -x * np.exp(-0.5 * x*x)
+        return 1./np.sqrt(2*np.pi) * -x * np.exp(-0.5 * x*x)
 
     def cdf(self, x):
         return special.ndtr(x)
@@ -181,6 +176,33 @@ def time_pinv_rvs(self, dist):
         self.rng.rvs(100000)
 
 
+class NumericalInverseHermite(Benchmark):
+
+    param_names = ['dist', 'order']
+    params = [allcontdists, [3, 5]]
+
+    def setup(self, dist, order):
+        self.urng = np.random.default_rng(0xb235b58c1f616c59c18d8568f77d44d1)
+        with np.testing.suppress_warnings() as sup:
+            sup.filter(RuntimeWarning)
+            try:
+                self.rng = stats.NumericalInverseHermite(
+                    dist, order=order, random_state=self.urng
+                )
+            except stats.UNURANError:
+                raise NotImplementedError(f"setup failed for {dist}")
+
+    def time_hinv_setup(self, dist, order):
+        with np.testing.suppress_warnings() as sup:
+            sup.filter(RuntimeWarning)
+            stats.NumericalInverseHermite(
+                dist, order=order, random_state=self.urng
+            )
+
+    def time_hinv_rvs(self, dist, order):
+        self.rng.rvs(100000)
+
+
 class DiscreteAliasUrn(Benchmark):
 
     param_names = ['distribution']
diff --git a/doc/source/tutorial/stats/sampling.rst b/doc/source/tutorial/stats/sampling.rst
index 5ca53df3bda0..057dc95c36f3 100644
--- a/doc/source/tutorial/stats/sampling.rst
+++ b/doc/source/tutorial/stats/sampling.rst
@@ -83,7 +83,7 @@ different methods is shown in the table below.
 Methods for continuous distributions   Required Inputs  Optional Inputs  Setup Speed  Sampling Speed
 =====================================  ===============  ===============  ===========  ==============
 :class:`~TransformedDensityRejection`  pdf, dpdf        None             slow         fast
-:class:`~NumericalInverseHermite`      cdf              None             (very) slow  (very) fast
+:class:`~NumericalInverseHermite`      cdf              pdf, dpdf        (very) slow  (very) fast
 :class:`~NumericalInversePolynomial`   pdf              cdf              (very) slow  (very) fast
 =====================================  ===============  ===============  ===========  ==============
 
@@ -286,6 +286,7 @@ Generators in :mod:`scipy.stats`
    sampling_tdr
    sampling_dau
    sampling_pinv
+   sampling_hinv
 
 
 References
diff --git a/doc/source/tutorial/stats/sampling_hinv.rst b/doc/source/tutorial/stats/sampling_hinv.rst
new file mode 100644
index 000000000000..c4f830f8ccd9
--- /dev/null
+++ b/doc/source/tutorial/stats/sampling_hinv.rst
@@ -0,0 +1,191 @@
+.. _sampling-hinv:
+
+Hermite interpolation based INVersion of CDF (HINV)
+===================================================
+
+.. currentmodule:: scipy.stats
+
+* Required: CDF
+* Optional: PDF, dPDF
+* Speed:
+
+  * Set-up: (very) slow
+  * Sampling: (very) fast
+
+HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+in each interval the inverse CDF is approximated by polynomials constructed by means
+of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+the accuracy by splitting a particular interval without recomputations in unaffected
+intervals. Three types of splines are implemented: linear, cubic, and quintic
+interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+also requires PDF and quintic interpolation PDF and its derivative.
+
+These splines have to be computed in a setup step. However, it only works for
+distributions with bounded domain; for distributions with unbounded domain the tails
+are chopped off such that the probability for the tail regions is small compared to
+the given u-resolution.
+
+The method is not exact, as it only produces random variates of the approximated
+distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
+quantile ``U`` i.e. ``X = approx_ppf(U)``) can be set to the
+required resolution (within machine precision). Notice that very small values of
+the u-resolution are possible but may increase the cost for the setup step.
+
+`NumericalInverseHermite` approximates the inverse of a continuous
+statistical distribution's CDF with a Hermite spline. Order of the
+hermite spline can be specified by passing the ``order`` parameter.
+
+As described in [1]_, it begins by evaluating the distribution's PDF and
+CDF at a mesh of quantiles ``x`` within the distribution's support.
+It uses the results to fit a Hermite spline ``H`` such that
+``H(p) == x``, where ``p`` is the array of percentiles corresponding
+with the quantiles ``x``. Therefore, the spline approximates the inverse
+of the distribution's CDF to machine precision at the percentiles ``p``,
+but typically, the spline will not be as accurate at the midpoints between
+the percentile points::
+
+    p_mid = (p[:-1] + p[1:])/2
+
+so the mesh of quantiles is refined as needed to reduce the maximum
+"u-error"::
+
+    u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
+
+below the specified tolerance `u_resolution`. Refinement stops when the required
+tolerance is achieved or when the number of mesh intervals after the next
+refinement could exceed the maximum allowed number `max_intervals`.
+
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+
+To create a generator to sample from the standard normal distribution, do:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+
+The `NumericalInverseHermite` has a method that approximates the PPF of the
+distribution.
+
+    >>> rng = NumericalInverseHermite(dist)
+    >>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
+    >>> np.allclose(rng.ppf(p), norm.ppf(p))
+    True
+
+Depending on the implementation of the distribution's random sampling
+method, the random variates generated may be nearly identical, given
+the same random state.
+
+    >>> dist = genexpon(9, 16, 3)
+    >>> rng = NumericalInverseHermite(dist)
+    >>> # `seed` ensures identical random streams are used by each `rvs` method
+    >>> seed = 500072020
+    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> np.allclose(rvs1, rvs2)
+    True
+
+To check that the random variates closely follow the given distribution, we can
+look at its histogram:
+
+.. plot::
+
+    >>> import numpy as np
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+    >>> import matplotlib.pyplot as plt
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+    >>> rvs = rng.rvs(10000)
+    >>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
+    >>> fx = norm.pdf(x)
+    >>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
+    >>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
+    >>> plt.xlabel('x')
+    >>> plt.ylabel('PDF(x)')
+    >>> plt.title('Numerical Inverse Hermite Samples')
+    >>> plt.legend()
+    >>> plt.show()
+
+Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use
+quintic Hermite interpolation to approximate the PPF by passing the ``order``
+parameter:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def dpdf(self, x):
+    ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
+
+Higher orders result is fewer number of intervals:
+
+    >>> rng3 = NumericalInverseHermite(dist, order=3)
+    >>> rng5 = NumericalInverseHermite(dist, order=5)
+    >>> rng3.intervals, rng5.intervals
+    (3000, 522)
+
+The u-error can be estimated by calling the `u_error` method. It runs a small
+Monte-Carlo simulation to estimate the u-error. By default, 100000 samples are
+used. To change this, you can pass the number of samples as an argument:
+
+    >>> rng1 = NumericalInverseHermite(dist, u_resolution=1.e-10)
+    >>> rng1.u_error(sample_size=1000000)  # uses one million samples
+    UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)
+
+This returns a namedtuple which contains the maximum u-error and the mean
+absolute u-error.
+
+The u-error can be reduced by decreasing the u-resolution (maximum allowed u-error):
+
+    >>> rng2 = NumericalInverseHermite(dist, u_resolution=1.e-13)
+    >>> rng2.u_error(sample_size=1000000)
+    UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)
+
+Note that this comes at the cost of increased setup time and number of intervals.
+
+    >>> rng1.intervals
+    1022
+    >>> rng2.intervals
+    5687
+    >>> from timeit import timeit
+    >>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-10)
+    >>> timeit(f, number=1)
+    0.017409582000254886  # may vary
+    >>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)
+    >>> timeit(f, number=1)
+    0.08671202100003939  # may vary
+
+Please see [1]_ and [2]_ for more details on this method.
+
+References
+----------
+
+.. [1] Hörmann, Wolfgang, and Josef Leydold. "Continuous random variate
+       generation by fast numerical inversion." ACM Transactions on
+       Modeling and Computer Simulation (TOMACS) 13.4 (2003): 347-362.
+.. [2] UNU.RAN reference manual, Section 5.3.5,
+       "HINV - Hermite interpolation based INVersion of CDF",
+       https://statmath.wu.ac.at/software/unuran/doc/unuran.html#HINV
diff --git a/doc/source/tutorial/stats/sampling_pinv.rst b/doc/source/tutorial/stats/sampling_pinv.rst
index a152f25bd915..1b63662dddba 100644
--- a/doc/source/tutorial/stats/sampling_pinv.rst
+++ b/doc/source/tutorial/stats/sampling_pinv.rst
@@ -209,7 +209,7 @@ more intervals compared to `NumericalInversePolynomial`:
     >>> # u-resolution of the generator.
     >>> rng_hermite = NumericalInverseHermite(norm(), tol=1e-12)
     >>> rng_hermite.intervals
-    3340
+    3000
     >>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
     >>> rng_poly.intervals
     252
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index 51b7d6acbba2..e86169c13aaa 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -456,7 +456,7 @@
 from ._bootstrap import bootstrap
 from ._entropy import *
 from ._hypotests import *
-from ._rvs_sampling import rvs_ratio_uniforms, NumericalInverseHermite
+from ._rvs_sampling import rvs_ratio_uniforms  # noqa
 from ._unuran import *  # noqa
 from ._page_trend_test import page_trend_test
 from ._mannwhitneyu import mannwhitneyu
diff --git a/scipy/stats/_rvs_sampling.py b/scipy/stats/_rvs_sampling.py
index 39a204e27a24..7df371006dbc 100644
--- a/scipy/stats/_rvs_sampling.py
+++ b/scipy/stats/_rvs_sampling.py
@@ -1,8 +1,6 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 from scipy._lib._util import check_random_state
-from scipy.interpolate import CubicHermiteSpline
-from scipy import stats
 
 
 def rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=1, c=0, random_state=None):
@@ -171,451 +169,3 @@ def rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=1, c=0, random_state=None):
         i += 1
 
     return np.reshape(x, size1d)
-
-
-class NumericalInverseHermite:
-    r"""
-    A Hermite spline fast numerical inverse of a probability distribution.
-
-    The initializer of `NumericalInverseHermite` accepts `dist`, an object
-    representing a continuous distribution, and provides an object with methods
-    that approximate `dist.ppf` and `dist.rvs`. For most distributions,
-    these methods are faster than those of `dist` itself.
-
-    Parameters
-    ----------
-    dist : object
-        Object representing the distribution for which a fast numerical inverse
-        is desired; for instance, a frozen instance of a `scipy.stats`
-        continuous distribution. See Notes and Examples for details.
-    tol : float, optional
-        u-error tolerance (see Notes). The default is 1e-12.
-    max_intervals : int, optional
-        Maximum number of intervals in the cubic Hermite spline used to
-        approximate the percent point function. The default is 100000.
-
-    Attributes
-    ----------
-    intervals : int
-        The number of intervals of the interpolant.
-    midpoint_error : float
-        The maximum u-error at an interpolant interval midpoint.
-
-    Notes
-    -----
-    `NumericalInverseHermite` approximates the inverse of a continuous
-    statistical distribution's CDF with a cubic Hermite spline.
-
-    As described in [1]_, it begins by evaluating the distribution's PDF and
-    CDF at a mesh of quantiles ``x`` within the distribution's support.
-    It uses the results to fit a cubic Hermite spline ``H`` such that
-    ``H(p) == x``, where ``p`` is the array of percentiles corresponding
-    with the quantiles ``x``. Therefore, the spline approximates the inverse
-    of the distribution's CDF to machine precision at the percentiles ``p``,
-    but typically, the spline will not be as accurate at the midpoints between
-    the percentile points::
-
-        p_mid = (p[:-1] + p[1:])/2
-
-    so the mesh of quantiles is refined as needed to reduce the maximum
-    "u-error"::
-
-        u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
-
-    below the specified tolerance `tol`. Refinement stops when the required
-    tolerance is achieved or when the number of mesh intervals after the next
-    refinement could exceed the maximum allowed number `max_intervals`.
-
-    The object `dist` must have methods ``pdf``, ``cdf``, and ``ppf`` that
-    behave like those of a *frozen* instance of `scipy.stats.rv_continuous`.
-    Specifically, it must have methods ``pdf`` and ``cdf`` that accept exactly
-    one ndarray argument ``x`` and return the probability density function and
-    cumulative density function (respectively) at ``x``. The object must also
-    have a method ``ppf`` that accepts a float ``p`` and returns the percentile
-    point function at ``p``. The object may also have a method ``isf`` that
-    accepts a float ``p`` and returns the inverse survival function at ``p``;
-    if it does not, it will be assigned an attribute ``isf`` that calculates
-    the inverse survival function using ``ppf``. The ``ppf`` and
-    ``isf` methods will each be evaluated at a small positive float ``p``
-    (e.g. ``p = utol/10``), and the domain over which the approximate numerical
-    inverse is defined will be ``ppf(p)`` to ``isf(p)``. The approximation will
-    not be accurate in the extreme tails beyond this domain.
-
-    References
-    ----------
-    .. [1] Hörmann, Wolfgang, and Josef Leydold. "Continuous random variate
-           generation by fast numerical inversion." ACM Transactions on
-           Modeling and Computer Simulation (TOMACS) 13.4 (2003): 347-362.
-
-    Examples
-    --------
-    For some distributions, ``dist.ppf`` and ``dist.rvs`` are quite slow.
-    For instance, consider `scipy.stats.genexpon`. We freeze the distribution
-    by passing all shape parameters into its initializer and time the resulting
-    object's ``ppf`` and ``rvs`` functions.
-
-    >>> import numpy as np
-    >>> from scipy import stats
-    >>> from timeit import timeit
-    >>> time_once = lambda f: f"{timeit(f, number=1)*1000:.6} ms"
-    >>> dist = stats.genexpon(9, 16, 3)  # freeze the distribution
-    >>> p = np.linspace(0.01, 0.99, 99)  # percentiles from 1% to 99%
-    >>> time_once(lambda: dist.ppf(p))
-    '154.565 ms'  # may vary
-
-    >>> time_once(lambda: dist.rvs(size=100))
-    '148.979 ms'  # may vary
-
-    The `NumericalInverseHermite` has a method that approximates ``dist.ppf``.
-
-    >>> from scipy.stats import NumericalInverseHermite
-    >>> fni = NumericalInverseHermite(dist)
-    >>> np.allclose(fni.ppf(p), dist.ppf(p))
-    True
-
-    In some cases, it is faster to both generate the fast numerical inverse
-    and use it than to call ``dist.ppf``.
-
-    >>> def time_me():
-    ...     fni = NumericalInverseHermite(dist)
-    ...     fni.ppf(p)
-    >>> time_once(time_me)
-    '11.9222 ms'  # may vary
-
-    After generating the fast numerical inverse, subsequent calls to its
-    methods are much faster.
-    >>> time_once(lambda: fni.ppf(p))
-    '0.0819 ms'  # may vary
-
-    The fast numerical inverse can also be used to generate random variates
-    using inverse transform sampling.
-
-    >>> time_once(lambda: fni.rvs(size=100))
-    '0.0911 ms'  # may vary
-
-    Depending on the implementation of the distribution's random sampling
-    method, the random variates generated may be nearly identical, given
-    the same random state.
-
-    >>> # `seed` ensures identical random streams are used by each `rvs` method
-    >>> seed = 500072020
-    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
-    >>> rvs2 = fni.rvs(size=100, random_state=np.random.default_rng(seed))
-    >>> np.allclose(rvs1, rvs2)
-    True
-
-    To use `NumericalInverseHermite` with a custom distribution, users may
-    subclass  `scipy.stats.rv_continuous` and initialize a frozen instance or
-    create an object with equivalent ``pdf``, ``cdf``, and ``ppf`` methods.
-    For instance, the following object represents the standard normal
-    distribution. For simplicity, we use `scipy.special.ndtr` and
-    `scipy.special.ndtri` to compute the ``cdf`` and ``ppf``, respectively.
-
-    >>> from scipy.special import ndtr, ndtri
-    >>>
-    >>> class MyNormal:
-    ...
-    ...     def pdf(self, x):
-    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
-    ...
-    ...     def cdf(self, x):
-    ...        return ndtr(x)
-    ...
-    ...     def ppf(self, x):
-    ...        return ndtri(x)
-    ...
-    >>> dist1 = MyNormal()
-    >>> fni1 = NumericalInverseHermite(dist1)
-    >>>
-    >>> dist2 = stats.norm()
-    >>> fni2 = NumericalInverseHermite(dist2)
-    >>>
-    >>> print(fni1.rvs(random_state=seed), fni2.rvs(random_state=seed))
-    -1.9603810921759424 -1.9603810921747074
-
-    """
-
-    def __init__(self, dist, *, tol=1e-12, max_intervals=100000):
-        res = _fast_numerical_inverse(dist, tol, max_intervals)
-        H, eu, intervals, a, b = res
-        self.H = H
-        self.midpoint_error = eu
-        self.intervals = intervals
-        self._a, self._b = a, b
-
-    def ppf(self, q):
-        r"""
-        Approximate percent point function (inverse `cdf`) of the given RV.
-
-        Parameters
-        ----------
-        q : array_like
-            lower tail probability.
-
-        Returns
-        -------
-        x : array_like
-            quantile corresponding to the lower tail probability `q`.
-
-        """
-        q = np.asarray(q)  # no harm; self.H always returns an array
-        result = np.zeros_like(q, dtype=np.float64)
-        i = (q >= 0) & (q <= 1)
-        result[i] = self.H(q[i])
-        result[~i] = np.nan
-        return result
-
-    def rvs(self, size=None, random_state=None):
-        """
-        Random variates of the given RV.
-
-        The `random_state` is used to draw uniform pseudo-random variates, and
-        these are converted to pseudo-random variates of the given RV using
-        inverse transform sampling.
-
-        Parameters
-        ----------
-        size : int, tuple of ints, or None; optional
-            Defines shape of array of random variates. Default is ``None``.
-        random_state : {None, int, `numpy.random.Generator`,
-                        `numpy.random.RandomState`}, optional
-
-            Defines the object to use for drawing pseudorandom variates.
-            If `random_state` is ``None`` the `np.random.RandomState`
-            singleton is used.
-            If `random_state` is an ``int``, a new ``RandomState`` instance is
-            used, seeded with `random_state`.
-            If `random_state` is already a ``RandomState`` or ``Generator``
-            instance, then that object is used.
-            Default is None.
-
-        Returns
-        -------
-        rvs : ndarray or scalar
-            Random variates of given `size`. If `size` is ``None``, a scalar
-            is returned.
-        """
-        random_state = check_random_state(random_state)
-        uniform = random_state.uniform(size=size)
-        # scale to valid domain of interpolant
-        uniform = self._a + uniform * (self._b - self._a)
-        return self.ppf(uniform)
-
-    def qrvs(self, size=None, d=None, qmc_engine=None):
-        """
-        Quasi-random variates of the given RV.
-
-        The `qmc_engine` is used to draw uniform quasi-random variates, and
-        these are converted to quasi-random variates of the given RV using
-        inverse transform sampling.
-
-        Parameters
-        ----------
-        size : int, tuple of ints, or None; optional
-            Defines shape of random variates array. Default is ``None``.
-        d : int or None, optional
-            Defines dimension of uniform quasi-random variates to be
-            transformed. Default is ``None``.
-        qmc_engine : scipy.stats.qmc.QMCEngine(d=1), optional
-            Defines the object to use for drawing
-            quasi-random variates. Default is ``None``, which uses
-            `scipy.stats.qmc.Halton(1)`.
-
-        Returns
-        -------
-        rvs : ndarray or scalar
-            Quasi-random variates. See Notes for shape information.
-
-        Notes
-        -----
-        The shape of the output array depends on `size`, `d`, and `qmc_engine`.
-        The intent is for the interface to be natural, but the detailed rules
-        to achieve this are complicated.
-
-        - If `qmc_engine` is ``None``, a `scipy.stats.qmc.Halton` instance is
-          created with dimension `d`. If `d` is not provided, ``d=1``.
-        - If `qmc_engine` is not ``None`` and `d` is ``None``, `d` is
-          determined from the dimension of the `qmc_engine`.
-        - If `qmc_engine` is not ``None`` and `d` is not ``None`` but the
-          dimensions are inconsistent, a ``ValueError`` is raised.
-        - After `d` is determined according to the rules above, the output
-          shape is ``tuple_shape + d_shape``, where:
-
-              - ``tuple_shape = tuple()`` if `size` is ``None``,
-              - ``tuple_shape = (size,)`` if `size` is an ``int``,
-              - ``tuple_shape = size`` if `size` is a sequence,
-              - ``d_shape = tuple()`` if `d` is ``None`` or `d` is 1, and
-              - ``d_shape = (d,)`` if `d` is greater than 1.
-
-        The elements of the returned array are part of a low-discrepancy
-        sequence. If `d` is 1, this means that none of the samples are truly
-        independent. If `d` > 1, each slice ``rvs[..., i]`` will be of a
-        quasi-independent sequence; see `scipy.stats.qmc.QMCEngine` for
-        details. Note that when `d` > 1, the samples returned are still those
-        of the provided univariate distribution, not a multivariate
-        generalization of that distribution.
-
-        """
-        # Input validation for `qmc_engine` and `d`
-        # Error messages for invalid `d` are raised by QMCEngine
-        # we could probably use a stats.qmc.check_qrandom_state
-        if isinstance(qmc_engine, stats.qmc.QMCEngine):
-            message = "`d` must be consistent with dimension of `qmc_engine`."
-            if d is not None and qmc_engine.d != d:
-                raise ValueError(message)
-            d = qmc_engine.d if d is None else d
-        elif qmc_engine is None:
-            d = 1 if d is None else d
-            qmc_engine = stats.qmc.Halton(d)
-        else:
-            message = ("`qmc_engine` must be an instance of "
-                       "`scipy.stats.qmc.QMCEngine` or `None`.")
-            raise ValueError(message)
-
-        # `rvs` is flexible about whether `size` is an int or tuple, so this
-        # should be, too.
-        try:
-            tuple_size = tuple(size)
-        except TypeError:
-            tuple_size = (size,)
-
-        # Get uniform QRVS from qmc_random and transform it
-        uniform = qmc_engine.random(np.prod(tuple_size) or 1)
-        # scale to valid domain of interpolant
-        uniform = self._a + uniform * (self._b - self._a)
-        qrvs = self.ppf(uniform)
-
-        # Output reshaping for user convenience
-        if size is None:
-            return qrvs.squeeze()[()]
-        else:
-            if d == 1:
-                return qrvs.reshape(tuple_size)
-            else:
-                return qrvs.reshape(tuple_size + (d,))
-
-
-def _fni_input_validation(dist, tol, max_intervals):
-    """
-    Input validation and standardization for _fast_numerical_inverse.
-
-    """
-
-    has_pdf = hasattr(dist, 'pdf') and callable(dist.pdf)
-    has_cdf = hasattr(dist, 'cdf') and callable(dist.cdf)
-    has_ppf = hasattr(dist, 'ppf') and callable(dist.ppf)
-    has_isf = hasattr(dist, 'isf') and callable(dist.isf)
-
-    if not (has_pdf and has_cdf and has_ppf):
-        raise ValueError("`dist` must have methods `pdf`, `cdf`, and `ppf`.")
-
-    if not has_isf:
-        def isf(x):
-            return 1 - dist.ppf(x)
-        dist.isf = isf
-
-    tol = float(tol)  # if there's an exception, raise it now
-
-    if int(max_intervals) != max_intervals or max_intervals <= 1:
-        raise ValueError("`max_intervals' must be an integer greater than 1.")
-
-    return dist, tol, max_intervals
-
-
-def _fast_numerical_inverse(dist, tol=1e-12, max_intervals=100000):
-    """
-    Generate fast, approximate PPF (inverse CDF) of probability distribution.
-
-    `_fast_numerical_inverse` accepts `dist`, an object representing the
-    distribution for which a fast approximate PPF is desired, and returns an
-    object `fni` with methods that approximate `dist.ppf` and `dist.rvs`.
-    For some distributions, these methods may be faster than those of `dist`
-    itself.
-
-    Parameters
-    ----------
-    dist : object
-        Object representing distribution for which fast approximate PPF is
-        desired; e.g., a frozen instance of `scipy.stats.rv_continuous`.
-    tol : float, optional
-        u-error tolerance. The default is 1e-12.
-    max_intervals : int, optional
-        Maximum number of intervals in the cubic Hermite Spline used to
-        approximate the percent point function. The default is 100000.
-
-    Returns
-    -------
-    H : scipy.interpolate.CubicHermiteSpline
-        Interpolant of the distributions's PPF.
-    intervals : int
-        The number of intervals of the interpolant.
-    midpoint_error : float
-        The maximum u-error at an interpolant interval midpoint.
-    a, b : float
-        The left and right endpoints of the valid domain of the interpolant.
-
-    """
-    dist, tol, max_intervals = _fni_input_validation(dist, tol, max_intervals)
-
-    # [1] Section 2.1: "For distributions with unbounded domain, we have to
-    # chop off its tails at [a] and [b] such that F(a) and 1-F(b) are small
-    # compared to the maximal tolerated approximation error."
-    p = np.array([dist.ppf(tol/10), dist.isf(tol/10)])  # initial interval
-
-    # [1] Section 2.3: "We then halve this interval recursively until
-    # |u[i+1]-u[i]| is smaller than some threshold value, for example, 0.05."
-    u = dist.cdf(p)
-    while p.size-1 <= np.ceil(max_intervals/2):
-        i = np.nonzero(np.diff(u) > 0.05)[0]
-        if not i.size:
-            break
-
-        p_mid = (p[i] + p[i+1])/2
-        # Compute only the new values and insert them in the right places
-        # [1] uses a linked list; we can't do that efficiently
-        u_mid = dist.cdf(p_mid)
-        p = np.concatenate((p, p_mid))
-        u = np.concatenate((u, u_mid))
-        i_sort = np.argsort(p)
-        p = p[i_sort]
-        u = u[i_sort]
-
-    # [1] Section 2.3: "Now we continue with checking the error estimate in
-    # each of the intervals and continue with splitting them until [it] is
-    # smaller than a given error bound."
-    u = dist.cdf(p)
-    f = dist.pdf(p)
-    while p.size-1 <= max_intervals:
-        # [1] Equation 4-8
-        try:
-            H = CubicHermiteSpline(u, p, 1/f)
-        except ValueError:
-            message = ("The interpolating spline could not be created. This "
-                       "is often caused by inaccurate CDF evaluation in a "
-                       "tail of the distribution. Increasing `tol` can "
-                       "resolve this error at the expense of lower accuracy.")
-            raise ValueError(message)
-        # To improve performance, add update feature to CubicHermiteSpline
-
-        # [1] Equation 12
-        u_mid = (u[:-1] + u[1:])/2
-        eu = np.abs(dist.cdf(H(u_mid)) - u_mid)
-
-        i = np.nonzero(eu > tol)[0]
-        if not i.size:
-            break
-
-        p_mid = (p[i] + p[i+1])/2
-        u_mid = dist.cdf(p_mid)
-        f_mid = dist.pdf(p_mid)
-        p = np.concatenate((p, p_mid))
-        u = np.concatenate((u, u_mid))
-        f = np.concatenate((f, f_mid))
-        i_sort = np.argsort(p)
-        p = p[i_sort]
-        u = u[i_sort]
-        f = f[i_sort]
-
-    # todo: add test for monotonicity [1] Section 2.4
-    # todo: deal with vanishing density [1] Section 2.5
-    return H, eu, p.size-1, u[0], u[-1]
diff --git a/scipy/stats/_unuran/__init__.py b/scipy/stats/_unuran/__init__.py
index 6974f2d02619..8638754a87eb 100644
--- a/scipy/stats/_unuran/__init__.py
+++ b/scipy/stats/_unuran/__init__.py
@@ -2,5 +2,6 @@
     TransformedDensityRejection,
     DiscreteAliasUrn,
     NumericalInversePolynomial,
+    NumericalInverseHermite,
     UNURANError
 )
diff --git a/scipy/stats/_unuran/unuran_wrapper.pyi b/scipy/stats/_unuran/unuran_wrapper.pyi
index 11743df2f13d..b5391dec9106 100644
--- a/scipy/stats/_unuran/unuran_wrapper.pyi
+++ b/scipy/stats/_unuran/unuran_wrapper.pyi
@@ -3,6 +3,7 @@ from typing import Union, Any, Tuple, List, overload, Callable, NamedTuple
 from typing_extensions import Protocol
 import numpy.typing as npt
 from scipy._lib._util import SeedType
+import scipy.stats as stats
 
 
 ArrayLike0D = Union[bool, int, float, complex, str, bytes, np.generic]
@@ -92,6 +93,39 @@ class NumericalInversePolynomial(Method):
     def u_error(self, sample_size: int = ...) -> UError: ...
 
 
+class HINVDist(Protocol):
+    @property
+    def pdf(self) -> Callable[..., float]: ...
+    @property
+    def cdf(self) -> Callable[..., float]: ...
+    @property
+    def support(self) -> Tuple[float, float]: ...
+
+
+class NumericalInverseHermite(Method):
+    def __init__(self,
+                 dist: HINVDist,
+                 *,
+                 domain: None | Tuple[float, float] = ...,
+                 order: int= ...,
+                 u_resolution: float = ...,
+                 construction_points: None | npt.ArrayLike = ...,
+                 max_intervals: int = ...,
+                 guide_factor: float = ...,
+                 random_state: SeedType = ...) -> None: ...
+    @property
+    def intervals(self) -> int: ...
+    @overload
+    def ppf(self, u: ArrayLike0D) -> float: ...  # type: ignore[misc]
+    @overload
+    def ppf(self, u: npt.ArrayLike) -> np.ndarray: ...
+    def qrvs(self,
+             size: None | int | Tuple[int, ...] = ...,
+             d: None | int = ...,
+             qmc_engine: None | stats.qmc.QMCEngine = ...) -> npt.ArrayLike: ...
+    def u_error(self, sample_size: int = ...) -> UError: ...
+
+
 class DAUDist(Protocol):
     @property
     def pmf(self) -> Callable[..., float]: ...
diff --git a/scipy/stats/_unuran/unuran_wrapper.pyx.templ b/scipy/stats/_unuran/unuran_wrapper.pyx.templ
index 52647751939a..441b2212d798 100644
--- a/scipy/stats/_unuran/unuran_wrapper.pyx.templ
+++ b/scipy/stats/_unuran/unuran_wrapper.pyx.templ
@@ -15,10 +15,12 @@ IF not NPY_OLD:
 from scipy._lib.ccallback cimport ccallback_t
 from scipy._lib.messagestream cimport MessageStream
 from .unuran cimport *
+import warnings
 import threading
 import functools
 from collections import namedtuple
 import numpy as np
+import scipy.stats as stats
 from scipy.stats._distn_infrastructure import argsreduce
 from scipy._lib._util import check_random_state
 
@@ -57,6 +59,8 @@ class UNURANError(RuntimeError):
     pass
 
 
+ctypedef double (*URNG_FUNCT)(void *) nogil
+
 IF not NPY_OLD:
     cdef object get_numpy_rng(object seed = None):
         """
@@ -78,8 +82,6 @@ IF not NPY_OLD:
             return np.random.default_rng(seed._bit_generator)
         return seed
 ELSE:
-    ctypedef double (*URNG_FUNCT)(void *) nogil
-
     cdef object get_numpy_rng(object seed = None):
         """
         Create a NumPy RandomState object from a given seed. If the seed is
@@ -112,6 +114,8 @@ cdef class _URNG:
         random number generator.
     """
     cdef object numpy_rng
+    cdef double[::1] qrvs_array
+    cdef size_t i
 
     def __init__(self, numpy_rng):
         self.numpy_rng = numpy_rng
@@ -121,6 +125,12 @@ cdef class _URNG:
             with gil:
                 return self.numpy_rng.uniform()
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef double _next_qdouble(self) nogil:
+        self.i += 1
+        return self.qrvs_array[self.i-1]
+
     cdef unur_urng * get_urng(self) except *:
         """
         Get a ``unur_urng`` object from given ``numpy_rng``.
@@ -151,12 +161,20 @@ cdef class _URNG:
 
             return unuran_urng
 
+    cdef unur_urng *get_qurng(self, size, qmc_engine) except *:
+        cdef unur_urng *unuran_urng
+        self.i = 0
+        self.qrvs_array = np.ascontiguousarray(qmc_engine.random(size).ravel())
+        unuran_urng = unur_urng_new(<URNG_FUNCT>self._next_qdouble,
+                                    <void *>self)
+        return unuran_urng
+
 
 # Module level lock. This is used to provide thread-safe error reporting.
 # UNU.RAN has a thread-unsafe global FILE streams where errors are logged.
 # To make it thread-safe, one can aquire a lock before calling
 # `unur_set_stream` and release once the stream is not needed anymore.
-cdef object _lock = threading.Lock()
+cdef object _lock = threading.RLock()
 
 cdef:
     unur_urng *default_urng
@@ -458,9 +476,9 @@ cdef class Method:
             if self._callback_wrapper is not None:
                 release_unuran_callback(&callback)
 
-    def rvs(self, size=None):
+    def rvs(self, size=None, random_state=None):
         """
-        rvs(size=None)
+        rvs(size=None, random_state=None)
 
         Sample from the distribution.
 
@@ -469,6 +487,17 @@ cdef class Method:
         size : int or tuple, optional
             The shape of samples. Default is ``None`` in which case a scalar
             sample is returned.
+        random_state : {None, int, `numpy.random.Generator`,
+                        `numpy.random.RandomState`}, optional
+
+            A NumPy random number generator or seed for the underlying NumPy random
+            number generator used to generate the stream of uniform random numbers.
+            If `random_state` is None (or `np.random`), `random_state` provided during
+            initialization is used.
+            If `random_state` is an int, a new ``RandomState`` instance is used,
+            seeded with `random_state`.
+            If `random_state` is already a ``Generator`` or ``RandomState`` instance then
+            that instance is used.
 
         Returns
         -------
@@ -478,15 +507,22 @@ cdef class Method:
         cdef double[::1] out_cont
         cdef int[::1] out_discr
         N = 1 if size is None else np.prod(size)
+        prev_random_state = self.numpy_rng
+        if random_state is not None:
+            self.set_random_state(random_state)
         if unur_distr_is_cont(unur_get_distr(self.rng)):
             out_cont = np.empty(N, dtype=np.float64)
             self._rvs_cont(out_cont)
+            if random_state is not None:
+                self.set_random_state(prev_random_state)
             if size is None:
                 return out_cont[0]
             return np.asarray(out_cont).reshape(size)
         elif unur_distr_is_discr(unur_get_distr(self.rng)):
             out_discr = np.empty(N, dtype=np.int32)
             self._rvs_discr(out_discr)
+            if random_state is not None:
+                self.set_random_state(prev_random_state)
             if size is None:
                 return out_discr[0]
             return np.asarray(out_discr).reshape(size)
@@ -1302,7 +1338,7 @@ cdef class NumericalInversePolynomial(Method):
         Parameters
         ----------
         u : array_like
-            Qunatiles.
+            Quantiles.
 
         Returns
         -------
@@ -1346,7 +1382,7 @@ cdef class NumericalInversePolynomial(Method):
         max_error : float
             Maximum u-error.
         mean_absolute_error : float
-            Maximum absolute u-error.
+            Mean absolute u-error.
         """
         # UNU.RAN doesn't return a proper error code for this condition.
         if sample_size < 1000:
@@ -1370,6 +1406,538 @@ cdef class NumericalInversePolynomial(Method):
         return UError(max_error, mae)
 
 
+cdef class NumericalInverseHermite(Method):
+    """
+    NumericalInverseHermite(dist, *, domain=None, order=3, u_resolution=1e-12, construction_points=None, max_intervals=100000, guide_factor=1, random_state=None)
+
+    Hermite interpolation based INVersion of CDF (HINV).
+
+    HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+    Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+    in each interval the inverse CDF is approximated by polynomials constructed by means
+    of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+    the accuracy by splitting a particular interval without recomputations in unaffected
+    intervals. Three types of splines are implemented: linear, cubic, and quintic
+    interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+    also requires PDF and quintic interpolation PDF and its derivative.
+
+    These splines have to be computed in a setup step. However, it only works for
+    distributions with bounded domain; for distributions with unbounded domain the tails
+    are chopped off such that the probability for the tail regions is small compared to
+    the given u-resolution.
+
+    The method is not exact, as it only produces random variates of the approximated
+    distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+    ``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
+    quantile ``U`` i.e. ``X = approx_ppf(U)``) can be set to the
+    required resolution (within machine precision). Notice that very small values of
+    the u-resolution are possible but may increase the cost for the setup step.
+
+    Parameters
+    ----------
+    dist : object
+        An instance of a class with a ``cdf`` and optionally a ``pdf`` and ``dpdf`` method.
+
+        * ``cdf``: CDF of the distribution. The signature of the CDF is expected to be:
+          ``def cdf(self, x: float) -> float``. i.e. the CDF should accept a Python
+          float and return a Python float.
+        * ``pdf``: PDF of the distribution. This method is optional when ``order=1``.
+          Must have the same signature as the PDF.
+        * ``dpdf``: Derivative of the PDF w.r.t the variate (i.e. ``x``). This method is
+          optional with ``order=1`` or ``order=3``. Must have the same signature as the CDF.
+
+    domain : list or tuple of length 2, optional
+        The support of the distribution.
+        Default is ``None``. When ``None``:
+
+        * If a ``support`` method is provided by the distribution object
+          `dist`, it is used to set the domain of the distribution.
+        * Otherwise the support is assumed to be :math:`(-\infty, \infty)`.
+
+    order : int, default: ``3``
+        Set order of Hermite interpolation. Valid orders are 1, 3, and 5.
+        Valid orders are 1, 3, and 5. Notice that order greater than 1 requires the density
+        of the distribution, and order greater than 3 even requires the derivative of the
+        density. Using order 1 results for most distributions in a huge number of intervals
+        and is therefore not recommended. If the maximal error in u-direction is very small
+        (say smaller than 1.e-10), order 5 is recommended as it leads to considerably fewer
+        design points, as long there are no poles or heavy tails.
+    u_resolution : float, default: ``1e-12``
+        Set maximal tolerated u-error. Notice that the resolution of most uniform random
+        number sources is 2-32= 2.3e-10. Thus a value of 1.e-10 leads to an inversion
+        algorithm that could be called exact. For most simulations slightly bigger values
+        for the maximal error are enough as well. Default is 1e-12.
+    construction_points : array_like, optional
+        Set starting construction points (nodes) for Hermite interpolation. As the possible
+        maximal error is only estimated in the setup it may be necessary to set some
+        special design points for computing the Hermite interpolation to guarantee that the
+        maximal u-error can not be bigger than desired. Such points are points where the
+        density is not differentiable or has a local extremum.
+    max_intervals : int, default: ``100000``
+        Maximum number of intervals in the Hermite spline used to
+        approximate the percent point function. The default is 100000.
+    guide_factor : float, default: ``1``
+        Set factor for relative size of the guide table for indexed search.
+        It must be greater than or equal to 0. When set to 0, then sequential
+        search is used. Default is 1.
+    random_state : {None, int, `numpy.random.Generator`,
+                        `numpy.random.RandomState`}, optional
+
+        A NumPy random number generator or seed for the underlying NumPy random
+        number generator used to generate the stream of uniform random numbers.
+        If `random_state` is None (or `np.random`), the `numpy.random.RandomState`
+        singleton is used.
+        If `random_state` is an int, a new ``RandomState`` instance is used,
+        seeded with `random_state`.
+        If `random_state` is already a ``Generator`` or ``RandomState`` instance then
+        that instance is used.
+
+    Notes
+    -----
+    `NumericalInverseHermite` approximates the inverse of a continuous
+    statistical distribution's CDF with a Hermite spline. Order of the
+    hermite spline can be specified by passing the ``order`` parameter.
+
+    As described in [1]_, it begins by evaluating the distribution's PDF and
+    CDF at a mesh of quantiles ``x`` within the distribution's support.
+    It uses the results to fit a Hermite spline ``H`` such that
+    ``H(p) == x``, where ``p`` is the array of percentiles corresponding
+    with the quantiles ``x``. Therefore, the spline approximates the inverse
+    of the distribution's CDF to machine precision at the percentiles ``p``,
+    but typically, the spline will not be as accurate at the midpoints between
+    the percentile points::
+
+        p_mid = (p[:-1] + p[1:])/2
+
+    so the mesh of quantiles is refined as needed to reduce the maximum
+    "u-error"::
+
+        u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
+
+    below the specified tolerance `u_resolution`. Refinement stops when the required
+    tolerance is achieved or when the number of mesh intervals after the next
+    refinement could exceed the maximum allowed number `max_intervals`.
+
+    This method accepted a `tol` parameter to specify the maximum tolerable u-error
+    but it has now been doc-deprecated in the favor of the `u_resolution` parameter.
+    A deprecation warning will be raised in v1.9.0 and the parameter will be removed
+    in v1.10.0.
+
+    References
+    ----------
+    .. [1] Hörmann, Wolfgang, and Josef Leydold. "Continuous random variate
+           generation by fast numerical inversion." ACM Transactions on
+           Modeling and Computer Simulation (TOMACS) 13.4 (2003): 347-362.
+    .. [2] UNU.RAN reference manual, Section 5.3.5,
+           "HINV - Hermite interpolation based INVersion of CDF",
+           https://statmath.wu.ac.at/software/unuran/doc/unuran.html#HINV
+
+    Examples
+    --------
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+
+    To create a generator to sample from the standard normal distribution, do:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+
+    The `NumericalInverseHermite` has a method that approximates the PPF of the
+    distribution.
+
+    >>> rng = NumericalInverseHermite(dist)
+    >>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
+    >>> np.allclose(rng.ppf(p), norm.ppf(p))
+    True
+
+    Depending on the implementation of the distribution's random sampling
+    method, the random variates generated may be nearly identical, given
+    the same random state.
+
+    >>> dist = genexpon(9, 16, 3)
+    >>> rng = NumericalInverseHermite(dist)
+    >>> # `seed` ensures identical random streams are used by each `rvs` method
+    >>> seed = 500072020
+    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> np.allclose(rvs1, rvs2)
+    True
+
+    To check that the random variates closely follow the given distribution, we can
+    look at its histogram:
+
+    >>> import matplotlib.pyplot as plt
+    >>> dist = StandardNormal()
+    >>> rng = NumericalInverseHermite(dist)
+    >>> rvs = rng.rvs(10000)
+    >>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
+    >>> fx = norm.pdf(x)
+    >>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
+    >>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
+    >>> plt.xlabel('x')
+    >>> plt.ylabel('PDF(x)')
+    >>> plt.title('Numerical Inverse Hermite Samples')
+    >>> plt.legend()
+    >>> plt.show()
+
+    Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use
+    quintic Hermite interpolation to approximate the PPF by passing the ``order``
+    parameter:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def dpdf(self, x):
+    ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
+
+    Higher orders result is fewer number of intervals:
+
+    >>> rng3 = NumericalInverseHermite(dist, order=3)
+    >>> rng5 = NumericalInverseHermite(dist, order=5)
+    >>> rng3.intervals, rng5.intervals
+    (3000, 522)
+
+    The u-error can be estimated by calling the `u_error` method. It runs a small
+    Monte-Carlo simulation to estimate the u-error. By default, 100000 samples are
+    used. To change this, you can pass the number of samples as an argument:
+
+    >>> rng1 = NumericalInverseHermite(dist, u_resolution=1.e-10)
+    >>> rng1.u_error(sample_size=1000000)  # uses one million samples
+    UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)
+
+    This returns a namedtuple which contains the maximum u-error and the mean
+    absolute u-error.
+
+    The u-error can be reduced by decreasing the u-resolution (maximum allowed u-error):
+
+    >>> rng2 = NumericalInverseHermite(dist, u_resolution=1.e-13)
+    >>> rng2.u_error(sample_size=1000000)
+    UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)
+
+    Note that this comes at the cost of increased setup time and number of intervals.
+
+    >>> rng1.intervals
+    1022
+    >>> rng2.intervals
+    5687
+    >>> from timeit import timeit
+    >>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-10)
+    >>> timeit(f, number=1)
+    0.017409582000254886  # may vary
+    >>> f = lambda: NumericalInverseHermite(dist, u_resolution=1e-13)
+    >>> timeit(f, number=1)
+    0.08671202100003939  # may vary
+
+    Since the PPF of the normal distribution is available as a special function, we
+    can also check the x-error, i.e. the difference between the approximated PPF and
+    exact PPF::
+
+    >>> import matplotlib.pyplot as plt
+    >>> u = np.linspace(0.01, 0.99, 1000)
+    >>> approxppf = rng.ppf(u)
+    >>> exactppf = norm.ppf(u)
+    >>> error = np.abs(exactppf - approxppf)
+    >>> plt.plot(u, error)
+    >>> plt.xlabel('u')
+    >>> plt.ylabel('error')
+    >>> plt.title('Error between exact and approximated PPF (x-error)')
+    >>> plt.show()
+
+    """
+    cdef double[::1] construction_points_array
+
+    def __cinit__(self,
+                  dist,
+                  *,
+                  domain=None,
+                  order=3,
+                  u_resolution=1e-12,
+                  tol=None,
+                  construction_points=None,
+                  max_intervals=100000,
+                  guide_factor=1,
+                  random_state=None):
+        domain, order, u_resolution = self._validate_args(dist, domain, order, max_intervals,
+                                                          u_resolution, tol, construction_points)
+
+        # save all the arguments for pickling support
+        self._kwargs = {
+            'dist': dist,
+            'domain': domain,
+            'order': order,
+            'u_resolution': u_resolution,
+            'tol': tol,
+            'max_intervals': max_intervals,
+            'guide_factor': guide_factor,
+            'construction_points': construction_points,
+            'random_state': random_state
+        }
+
+        cdef:
+            unur_distr *distr
+            unur_par *par
+
+        self.callbacks = _unpack_dist(dist, "cont", meths=["cdf"], optional_meths=["pdf", "dpdf"])
+        def _callback_wrapper(x, name):
+            return self.callbacks[name](x)
+        self._callback_wrapper = _callback_wrapper
+        self._messages = MessageStream()
+        _lock.acquire()
+        try:
+            unur_set_stream(self._messages.handle)
+
+            self.distr = unur_distr_cont_new()
+            if self.distr == NULL:
+                raise UNURANError(self._messages.get())
+            _pack_distr(self.distr, self.callbacks)
+
+            if domain is not None:
+                self._check_errorcode(unur_distr_cont_set_domain(self.distr, domain[0],
+                                                                 domain[1]))
+
+            self.par = unur_hinv_new(self.distr)
+            if self.par == NULL:
+                raise UNURANError(self._messages.get())
+
+            self._check_errorcode(unur_hinv_set_order(self.par, order))
+            self._check_errorcode(unur_hinv_set_u_resolution(self.par, u_resolution))
+            self._check_errorcode(unur_hinv_set_max_intervals(self.par, max_intervals))
+            self._check_errorcode(unur_hinv_set_cpoints(self.par, &self.construction_points_array[0],
+                                                        len(self.construction_points_array)))
+            self._check_errorcode(unur_hinv_set_guidefactor(self.par, guide_factor))
+
+            self._set_rng(random_state)
+        finally:
+            _lock.release()
+
+    def _validate_args(self, dist, domain, order, max_intervals, u_resolution, tol,
+                       construction_points):
+        domain = _validate_domain(domain, dist)
+        # UNU.RAN raises warning and sets a default value. Prefer an error instead.
+        if order not in {1, 3, 5}:
+            raise ValueError("`order` must be either 1, 3, or 5.")
+        if int(max_intervals) != max_intervals or max_intervals <= 1:
+            raise ValueError("`max_intervals' must be an integer greater than 1.")
+        if tol is not None:
+            u_resolution = tol
+        u_resolution = float(u_resolution)
+        self.construction_points_array = np.ascontiguousarray(construction_points,
+                                                              dtype=np.float64)
+        if len(self.construction_points_array) == 0:
+            raise ValueError("`construction_points` must be a non-empty array.")
+        return (domain, order, u_resolution)
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef inline void _ppf(self, const double *u, double *out, size_t N):
+        cdef:
+            size_t i
+        for i in range(N):
+            out[i] = unur_hinv_eval_approxinvcdf(self.rng, u[i])
+
+    def ppf(self, u):
+        """
+        ppf(u)
+
+        Approximated PPF of the given distribution.
+
+        Parameters
+        ----------
+        u : array_like
+            Quantiles.
+
+        Returns
+        -------
+        ppf : array_like
+            Percentiles corresponding to given quantiles `u`.
+        """
+        u = np.asarray(u, dtype='d')
+        oshape = u.shape
+        u = u.ravel()
+        # UNU.RAN fills in ends of the support when u < 0 or u > 1 while
+        # SciPy fills in nans. Prefer SciPy behaviour.
+        cond0 = 0 <= u
+        cond1 = u <= 1
+        cond2 = cond0 & cond1
+        goodu = argsreduce(cond2, u)[0]
+        out = np.empty_like(u)
+        cdef double[::1] u_view = np.ascontiguousarray(goodu)
+        cdef double[::1] goodout = np.empty_like(u_view)
+        if cond2.any():
+            self._ppf(&u_view[0], &goodout[0], len(goodu))
+        np.place(out, cond2, goodout)
+        np.place(out, ~cond2, np.nan)
+        return np.asarray(out).reshape(oshape)[()]
+
+    def u_error(self, sample_size=100000):
+        """
+        u_error(sample_size=100000)
+
+        Estimate the u-error of the approximation using Monte Carlo simulation.
+        This is only available if the generator was initialized with a `dist`
+        object containing the implementation of the exact CDF under `cdf` method.
+
+        Parameters
+        ----------
+        sample_size : int, optional
+            Number of samples to use for the estimation. It must be greater than
+            or equal to 1000.
+
+        Returns
+        -------
+        max_error : float
+            Maximum u-error.
+        mean_absolute_error : float
+            Mean absolute u-error.
+        """
+        # UNU.RAN doesn't return a proper error code for this condition.
+        if sample_size < 1000:
+            raise ValueError("`sample_size` must be greater than or equal to 1000.")
+        cdef double max_error, mae
+        cdef ccallback_t callback
+        _lock.acquire()
+        try:
+            self._messages.clear()
+            unur_set_stream(self._messages.handle)
+            init_unuran_callback(&callback, self._callback_wrapper)
+            self._check_errorcode(unur_hinv_estimate_error(self.rng, sample_size,
+                                                           &max_error, &mae))
+        finally:
+            _lock.release()
+            release_unuran_callback(&callback)
+        return UError(max_error, mae)
+
+    def qrvs(self, size=None, d=None, qmc_engine=None):
+        """
+        qrvs(size=None, d=None, qmc_engine=None)
+
+        Quasi-random variates of the given RV.
+
+        The `qmc_engine` is used to draw uniform quasi-random variates, and
+        these are converted to quasi-random variates of the given RV using
+        inverse transform sampling.
+
+        Parameters
+        ----------
+        size : int, tuple of ints, or None; optional
+            Defines shape of random variates array. Default is ``None``.
+        d : int or None, optional
+            Defines dimension of uniform quasi-random variates to be
+            transformed. Default is ``None``.
+        qmc_engine : scipy.stats.qmc.QMCEngine(d=1), optional
+            Defines the object to use for drawing
+            quasi-random variates. Default is ``None``, which uses
+            `scipy.stats.qmc.Halton(1)`.
+
+        Returns
+        -------
+        rvs : ndarray or scalar
+            Quasi-random variates. See Notes for shape information.
+
+        Notes
+        -----
+        The shape of the output array depends on `size`, `d`, and `qmc_engine`.
+        The intent is for the interface to be natural, but the detailed rules
+        to achieve this are complicated.
+
+        - If `qmc_engine` is ``None``, a `scipy.stats.qmc.Halton` instance is
+          created with dimension `d`. If `d` is not provided, ``d=1``.
+        - If `qmc_engine` is not ``None`` and `d` is ``None``, `d` is
+          determined from the dimension of the `qmc_engine`.
+        - If `qmc_engine` is not ``None`` and `d` is not ``None`` but the
+          dimensions are inconsistent, a ``ValueError`` is raised.
+        - After `d` is determined according to the rules above, the output
+          shape is ``tuple_shape + d_shape``, where:
+
+              - ``tuple_shape = tuple()`` if `size` is ``None``,
+              - ``tuple_shape = (size,)`` if `size` is an ``int``,
+              - ``tuple_shape = size`` if `size` is a sequence,
+              - ``d_shape = tuple()`` if `d` is ``None`` or `d` is 1, and
+              - ``d_shape = (d,)`` if `d` is greater than 1.
+
+        The elements of the returned array are part of a low-discrepancy
+        sequence. If `d` is 1, this means that none of the samples are truly
+        independent. If `d` > 1, each slice ``rvs[..., i]`` will be of a
+        quasi-independent sequence; see `scipy.stats.qmc.QMCEngine` for
+        details. Note that when `d` > 1, the samples returned are still those
+        of the provided univariate distribution, not a multivariate
+        generalization of that distribution.
+
+        """
+        # Input validation for `qmc_engine` and `d`
+        # Error messages for invalid `d` are raised by QMCEngine
+        # we could probably use a stats.qmc.check_qrandom_state
+        if isinstance(qmc_engine, stats.qmc.QMCEngine):
+            message = "`d` must be consistent with dimension of `qmc_engine`."
+            if d is not None and qmc_engine.d != d:
+                raise ValueError(message)
+            d = qmc_engine.d if d is None else d
+        elif qmc_engine is None:
+            d = 1 if d is None else d
+            qmc_engine = stats.qmc.Halton(d)
+        else:
+            message = ("`qmc_engine` must be an instance of "
+                       "`scipy.stats.qmc.QMCEngine` or `None`.")
+            raise ValueError(message)
+
+        # `rvs` is flexible about whether `size` is an int or tuple, so this
+        # should be, too.
+        try:
+            if size is None:
+                tuple_size = (1, )
+            else:
+                tuple_size = tuple(size)
+        except TypeError:
+            tuple_size = (size,)
+
+        cdef unur_urng *unuran_urng
+        cdef double[::1] qrvs_view
+        N = 1 if size is None else np.prod(size)
+        N = N*d
+        qrvs_view = np.empty(N, dtype=np.float64)
+        _lock.acquire()
+        try:
+            # the call below must be under a lock
+            unuran_urng = self._urng_builder.get_qurng(size=N, qmc_engine=qmc_engine)
+            unur_chg_urng(self.rng, unuran_urng)
+            self._rvs_cont(qrvs_view)
+            self.set_random_state(self.numpy_rng)
+            qrvs = np.asarray(qrvs_view).reshape(tuple_size + (d,))
+        finally:
+            _lock.release()
+
+        # Output reshaping for user convenience
+        if size is None:
+            return qrvs.squeeze()[()]
+        else:
+            if d == 1:
+                return qrvs.reshape(tuple_size)
+            else:
+                return qrvs.reshape(tuple_size + (d,))
+
+    @property
+    def intervals(self):
+        """
+        Get number of nodes (design points) used for Hermite interpolation in the
+        generator object. The number of intervals is the number of nodes minus 1.
+        """
+        return unur_hinv_get_n_intervals(self.rng)
+
+
 cdef class DiscreteAliasUrn(Method):
     r"""
     DiscreteAliasUrn(dist, *, domain=None, urn_factor=1, random_state=None)
diff --git a/scipy/stats/tests/test_sampling.py b/scipy/stats/tests/test_sampling.py
index b5b0a6fb93ef..8c7ec4ffaa7a 100644
--- a/scipy/stats/tests/test_sampling.py
+++ b/scipy/stats/tests/test_sampling.py
@@ -2,19 +2,22 @@
 import threading
 import pickle
 import pytest
+from copy import deepcopy
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal, suppress_warnings
 from numpy.lib import NumpyVersion
 from scipy.stats import (
     TransformedDensityRejection,
     DiscreteAliasUrn,
-    NumericalInversePolynomial
+    NumericalInversePolynomial,
+    NumericalInverseHermite
 )
 from scipy.stats import UNURANError
 from scipy import stats
 from scipy import special
 from scipy.stats import chisquare, cramervonmises
 from scipy.stats._distr_params import distdiscrete, distcont
+from scipy._lib._util import check_random_state
 
 
 # common test data: this data can be shared between all the tests.
@@ -23,10 +26,11 @@
 # Normal distribution shared between all the continuous methods
 class StandardNormal:
     def pdf(self, x):
-        return np.exp(-0.5 * x*x)
+        # normalization constant needed for NumericalInverseHermite
+        return 1./np.sqrt(2.*np.pi) * np.exp(-0.5 * x*x)
 
     def dpdf(self, x):
-        return -x * np.exp(-0.5 * x*x)
+        return 1./np.sqrt(2.*np.pi) * -x * np.exp(-0.5 * x*x)
 
     def cdf(self, x):
         return special.ndtr(x)
@@ -52,7 +56,8 @@ def support(self):
 all_methods = [
     ("TransformedDensityRejection", {"dist": StandardNormal()}),
     ("DiscreteAliasUrn", {"dist": [0.02, 0.18, 0.8]}),
-    ("NumericalInversePolynomial", {"dist": StandardNormal()})
+    ("NumericalInversePolynomial", {"dist": StandardNormal()}),
+    ("NumericalInverseHermite", {"dist": StandardNormal()})
 ]
 
 
@@ -835,3 +840,300 @@ def pdf(self, x):
         msg = r"Exact CDF required but not found."
         with pytest.raises(ValueError, match=msg):
             rng.u_error()
+
+
+class TestNumericalInverseHermite:
+    #         /  (1 +sin(2 Pi x))/2  if |x| <= 1
+    # f(x) = <
+    #         \  0        otherwise
+    # Taken from UNU.RAN test suite (from file t_hinv.c)
+    class dist0:
+        def pdf(self, x):
+            return 0.5*(1. + np.sin(2.*np.pi*x))
+
+        def dpdf(self, x):
+            return np.pi*np.cos(2.*np.pi*x)
+
+        def cdf(self, x):
+            return (1. + 2.*np.pi*(1 + x) - np.cos(2.*np.pi*x)) / (4.*np.pi)
+
+        def support(self):
+            return -1, 1
+
+    #         /  Max(sin(2 Pi x)),0)Pi/2  if -1 < x <0.5
+    # f(x) = <
+    #         \  0        otherwise
+    # Taken from UNU.RAN test suite (from file t_hinv.c)
+    class dist1:
+        def pdf(self, x):
+            if (x <= -0.5):
+                return np.sin((2. * np.pi) * x) * 0.5 * np.pi
+            if (x < 0.):
+                return 0.
+            if (x <= 0.5):
+                return np.sin((2. * np.pi) * x) * 0.5 * np.pi
+
+        def dpdf(self, x):
+            if (x <= -0.5):
+                return np.cos((2. * np.pi) * x) * np.pi * np.pi
+            if (x < 0.):
+                return 0.
+            if (x <= 0.5):
+                return np.cos((2. * np.pi) * x) * np.pi * np.pi
+
+        def cdf(self, x):
+            if (x <= -0.5):
+                return 0.25 * (1 - np.cos((2. * np.pi) * x))
+            if (x < 0.):
+                return 0.5
+            if (x <= 0.5):
+                return 0.75 - 0.25 * np.cos((2. * np.pi) * x)
+
+        def support(self):
+            return -1, 0.5
+
+    dists = [dist0(), dist1()]
+
+    # exact mean and variance of the distributions in the list dists
+    mv0 = [-1/(2*np.pi), 1/3 - 1/(4*np.pi*np.pi)]
+    mv1 = [-1/4, 3/8-1/(2*np.pi*np.pi) - 1/16]
+    mvs = [mv0, mv1]
+
+    @pytest.mark.parametrize("dist, mv_ex",
+                             zip(dists, mvs))
+    @pytest.mark.parametrize("order", [3, 5])
+    def test_basic(self, dist, mv_ex, order):
+        rng = NumericalInverseHermite(dist, order=order, random_state=42)
+        check_cont_samples(rng, dist, mv_ex)
+
+    # test domains with inf + nan in them. need to write a custom test for
+    # this because not all methods support infinite tails.
+    @pytest.mark.parametrize("domain, err, msg", inf_nan_domains)
+    def test_inf_nan_domains(self, domain, err, msg):
+        with pytest.raises(err, match=msg):
+            NumericalInverseHermite(StandardNormal(), domain=domain)
+
+    @pytest.mark.xslow
+    @pytest.mark.parametrize(("distname", "shapes"), distcont)
+    def test_basic_all_scipy_dists(self, distname, shapes):
+        slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
+        fail_dists = {'beta', 'gausshyper', 'geninvgauss', 'ncf', 'nct',
+                      'norminvgauss', 'genhyperbolic', 'studentized_range',
+                      'vonmises', 'kappa4'}
+
+        if distname in slow_dists:
+            pytest.skip("Distribution is too slow")
+        if distname in fail_dists:
+            # specific reasons documented in gh-13319
+            # https://github.com/scipy/scipy/pull/13319#discussion_r626188955
+            pytest.xfail("Fails - usually due to inaccurate CDF/PDF")
+
+        np.random.seed(0)
+
+        dist = getattr(stats, distname)(*shapes)
+
+        with np.testing.suppress_warnings() as sup:
+            sup.filter(RuntimeWarning, "overflow encountered")
+            sup.filter(RuntimeWarning, "divide by zero")
+            sup.filter(RuntimeWarning, "invalid value encountered")
+            fni = stats.NumericalInverseHermite(dist)
+
+        x = np.random.rand(10)
+        p_tol = np.max(np.abs(dist.ppf(x)-fni.ppf(x))/np.abs(dist.ppf(x)))
+        u_tol = np.max(np.abs(dist.cdf(fni.ppf(x)) - x))
+
+        assert p_tol < 1e-8
+        assert u_tol < 1e-12
+
+    def test_input_validation(self):
+        match = r"`order` must be either 1, 3, or 5."
+        with pytest.raises(ValueError, match=match):
+            NumericalInverseHermite(StandardNormal(), order=2)
+
+        match = "`cdf` required but not found"
+        with pytest.raises(ValueError, match=match):
+            stats.NumericalInverseHermite("norm")
+
+        match = "could not convert string to float"
+        with pytest.raises(ValueError, match=match):
+            stats.NumericalInverseHermite(StandardNormal(), tol='ekki')
+
+        match = "`max_intervals' must be..."
+        with pytest.raises(ValueError, match=match):
+            stats.NumericalInverseHermite(StandardNormal(), max_intervals=-1)
+
+        match = "`qmc_engine` must be an instance of..."
+        with pytest.raises(ValueError, match=match):
+            fni = stats.NumericalInverseHermite(StandardNormal())
+            fni.qrvs(qmc_engine=0)
+
+        if NumpyVersion(np.__version__) >= '1.18.0':
+            # issues with QMCEngines and old NumPy
+            fni = stats.NumericalInverseHermite(StandardNormal())
+
+            match = "`d` must be consistent with dimension of `qmc_engine`."
+            with pytest.raises(ValueError, match=match):
+                fni.qrvs(d=3, qmc_engine=stats.qmc.Halton(2))
+
+    rngs = [None, 0, np.random.RandomState(0)]
+    if NumpyVersion(np.__version__) >= '1.18.0':
+        rngs.append(np.random.default_rng(0))  # type: ignore
+    sizes = [(None, tuple()), (8, (8,)), ((4, 5, 6), (4, 5, 6))]
+
+    @pytest.mark.parametrize('rng', rngs)
+    @pytest.mark.parametrize('size_in, size_out', sizes)
+    def test_RVS(self, rng, size_in, size_out):
+        dist = StandardNormal()
+        fni = stats.NumericalInverseHermite(dist)
+
+        rng2 = deepcopy(rng)
+        rvs = fni.rvs(size=size_in, random_state=rng)
+        if size_in is not None:
+            assert(rvs.shape == size_out)
+
+        if rng2 is not None:
+            rng2 = check_random_state(rng2)
+            uniform = rng2.uniform(size=size_in)
+            rvs2 = stats.norm.ppf(uniform)
+            assert_allclose(rvs, rvs2)
+
+    if NumpyVersion(np.__version__) >= '1.18.0':
+        qrngs = [None, stats.qmc.Sobol(1, seed=0), stats.qmc.Halton(3, seed=0)]
+    else:
+        qrngs = []
+    # `size=None` should not add anything to the shape, `size=1` should
+    sizes = [(None, tuple()), (1, (1,)), (4, (4,)),
+             ((4,), (4,)), ((2, 4), (2, 4))]  # type: ignore
+    # Neither `d=None` nor `d=1` should add anything to the shape
+    ds = [(None, tuple()), (1, tuple()), (3, (3,))]
+
+    @pytest.mark.parametrize('qrng', qrngs)
+    @pytest.mark.parametrize('size_in, size_out', sizes)
+    @pytest.mark.parametrize('d_in, d_out', ds)
+    def test_QRVS(self, qrng, size_in, size_out, d_in, d_out):
+        dist = StandardNormal()
+        fni = stats.NumericalInverseHermite(dist)
+
+        # If d and qrng.d are inconsistent, an error is raised
+        if d_in is not None and qrng is not None and qrng.d != d_in:
+            match = "`d` must be consistent with dimension of `qmc_engine`."
+            with pytest.raises(ValueError, match=match):
+                fni.qrvs(size_in, d=d_in, qmc_engine=qrng)
+            return
+
+        # Sometimes d is really determined by qrng
+        if d_in is None and qrng is not None and qrng.d != 1:
+            d_out = (qrng.d,)
+
+        shape_expected = size_out + d_out
+
+        qrng2 = deepcopy(qrng)
+        qrvs = fni.qrvs(size=size_in, d=d_in, qmc_engine=qrng)
+        if size_in is not None:
+            assert(qrvs.shape == shape_expected)
+
+        if qrng2 is not None:
+            uniform = qrng2.random(np.prod(size_in) or 1)
+            qrvs2 = stats.norm.ppf(uniform).reshape(shape_expected)
+            assert_allclose(qrvs, qrvs2, atol=1e-12)
+
+    def test_QRVS_size_tuple(self):
+        # QMCEngine samples are always of shape (n, d). When `size` is a tuple,
+        # we set `n = prod(size)` in the call to qmc_engine.random, transform
+        # the sample, and reshape it to the final dimensions. When we reshape,
+        # we need to be careful, because the _columns_ of the sample returned
+        # by a QMCEngine are "independent"-ish, but the elements within the
+        # columns are not. We need to make sure that this doesn't get mixed up
+        # by reshaping: qrvs[..., i] should remain "independent"-ish of
+        # qrvs[..., i+1], but the elements within qrvs[..., i] should be
+        # transformed from the same low-discrepancy sequence.
+        if NumpyVersion(np.__version__) <= '1.18.0':
+            pytest.skip("QMC doesn't play well with old NumPy")
+
+        dist = StandardNormal()
+        fni = stats.NumericalInverseHermite(dist)
+
+        size = (3, 4)
+        d = 5
+        qrng = stats.qmc.Halton(d, seed=0)
+        qrng2 = stats.qmc.Halton(d, seed=0)
+
+        uniform = qrng2.random(np.prod(size))
+
+        qrvs = fni.qrvs(size=size, d=d, qmc_engine=qrng)
+        qrvs2 = stats.norm.ppf(uniform)
+
+        for i in range(d):
+            sample = qrvs[..., i]
+            sample2 = qrvs2[:, i].reshape(size)
+            assert_allclose(sample, sample2, atol=1e-12)
+
+    def test_inaccurate_CDF(self):
+        # CDF function with inaccurate tail cannot be inverted; see gh-13319
+        # https://github.com/scipy/scipy/pull/13319#discussion_r626188955
+        shapes = (2.3098496451481823, 0.6268795430096368)
+        match = ("98 : one or more intervals very short; possibly due to "
+                 "numerical problems with a pole or very flat tail")
+
+        # fails with default tol
+        with pytest.warns(RuntimeWarning, match=match):
+            stats.NumericalInverseHermite(stats.beta(*shapes))
+
+        # no error with coarser tol
+        stats.NumericalInverseHermite(stats.beta(*shapes), tol=1e-8)
+
+    def test_custom_distribution(self):
+        dist1 = StandardNormal()
+        fni1 = stats.NumericalInverseHermite(dist1)
+
+        dist2 = stats.norm()
+        fni2 = stats.NumericalInverseHermite(dist2)
+
+        assert_allclose(fni1.rvs(random_state=0), fni2.rvs(random_state=0))
+
+    u = [
+        # check the correctness of the PPF for equidistant points between
+        # 0.02 and 0.98.
+        np.linspace(0., 1., num=10000),
+        # test the PPF method for empty arrays
+        [], [[]],
+        # test if nans and infs return nan result.
+        [np.nan], [-np.inf, np.nan, np.inf],
+        # test if a scalar is returned for a scalar input.
+        0,
+        # test for arrays with nans, values greater than 1 and less than 0,
+        # and some valid values.
+        [[np.nan, 0.5, 0.1], [0.2, 0.4, np.inf], [-2, 3, 4]]
+    ]
+
+    @pytest.mark.parametrize("u", u)
+    def test_ppf(self, u):
+        dist = StandardNormal()
+        rng = NumericalInverseHermite(dist, u_resolution=1e-12)
+        # Older versions of NumPy throw RuntimeWarnings for comparisons
+        # with nan.
+        with suppress_warnings() as sup:
+            sup.filter(RuntimeWarning, "invalid value encountered in greater")
+            sup.filter(RuntimeWarning, "invalid value encountered in "
+                                       "greater_equal")
+            sup.filter(RuntimeWarning, "invalid value encountered in less")
+            sup.filter(RuntimeWarning, "invalid value encountered in "
+                                       "less_equal")
+            res = rng.ppf(u)
+            expected = stats.norm.ppf(u)
+        assert_allclose(res, expected, rtol=1e-9, atol=3e-10)
+        assert res.shape == expected.shape
+
+    def test_u_error(self):
+        dist = StandardNormal()
+        rng = NumericalInverseHermite(dist, u_resolution=1e-10)
+        max_error, mae = rng.u_error()
+        assert max_error < 1e-10
+        assert mae <= max_error
+        with suppress_warnings() as sup:
+            # ignore warning about u-resolution being too small.
+            sup.filter(RuntimeWarning)
+            rng = NumericalInverseHermite(dist, u_resolution=1e-14)
+        max_error, mae = rng.u_error()
+        assert max_error < 1e-14
+        assert mae <= max_error
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index 9cfd2a1439c1..bc3a77bf78fb 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -23,8 +23,6 @@
 from numpy import array, arange, float32, float64, power
 import numpy as np
 
-from scipy._lib._util import check_random_state
-from scipy import special
 import scipy.stats as stats
 import scipy.stats.mstats as mstats
 import scipy.stats.mstats_basic as mstats_basic
@@ -34,7 +32,6 @@
 from .common_tests import check_named_results
 from scipy.sparse.sputils import matrix
 from scipy.spatial.distance import cdist
-from scipy.stats._distr_params import distcont
 from numpy.lib import NumpyVersion
 from scipy.stats.stats import (_broadcast_concatenate,
                                AlexanderGovernConstantInputWarning)
@@ -7174,188 +7171,6 @@ def test_dist_perm(self):
         assert_approx_equal(pvalue_dist, 0.001, significant=1)
 
 
-class TestNumericalInverseHermite:
-    @pytest.mark.parametrize(("distname", "shapes"), distcont)
-    def test_basic(self, distname, shapes):
-        slow_dists = {'ksone', 'kstwo', 'levy_stable', 'skewnorm'}
-        fail_dists = {'beta', 'gausshyper', 'geninvgauss', 'ncf', 'nct',
-                      'norminvgauss', 'genhyperbolic', 'studentized_range'}
-
-        if distname in slow_dists:
-            pytest.skip("Distribution is too slow")
-        if distname in fail_dists:
-            # specific reasons documented in gh-13319
-            # https://github.com/scipy/scipy/pull/13319#discussion_r626188955
-            pytest.xfail("Fails - usually due to inaccurate CDF/PDF")
-
-        np.random.seed(0)
-
-        dist = getattr(stats, distname)(*shapes)
-
-        with np.testing.suppress_warnings() as sup:
-            sup.filter(RuntimeWarning, "overflow encountered")
-            sup.filter(RuntimeWarning, "divide by zero")
-            sup.filter(RuntimeWarning, "invalid value encountered")
-            fni = stats.NumericalInverseHermite(dist)
-
-        x = np.random.rand(10)
-        p_tol = np.max(np.abs(dist.ppf(x)-fni.ppf(x))/np.abs(dist.ppf(x)))
-        u_tol = np.max(np.abs(dist.cdf(fni.ppf(x)) - x))
-
-        assert p_tol < 1e-8
-        assert u_tol < 1e-12
-
-    def test_input_validation(self):
-        match = "`dist` must have methods `pdf`, `cdf`, and `ppf`"
-        with pytest.raises(ValueError, match=match):
-            stats.NumericalInverseHermite("norm")
-
-        match = "could not convert string to float"
-        with pytest.raises(ValueError, match=match):
-            stats.NumericalInverseHermite(stats.norm(), tol='ekki')
-
-        match = "`max_intervals' must be..."
-        with pytest.raises(ValueError, match=match):
-            stats.NumericalInverseHermite(stats.norm(), max_intervals=-1)
-
-        match = "`qmc_engine` must be an instance of..."
-        with pytest.raises(ValueError, match=match):
-            fni = stats.NumericalInverseHermite(stats.norm())
-            fni.qrvs(qmc_engine=0)
-
-        if NumpyVersion(np.__version__) >= '1.18.0':
-            # issues with QMCEngines and old NumPy
-            fni = stats.NumericalInverseHermite(stats.norm())
-
-            match = "`d` must be consistent with dimension of `qmc_engine`."
-            with pytest.raises(ValueError, match=match):
-                fni.qrvs(d=3, qmc_engine=stats.qmc.Halton(2))
-
-    rngs = [None, 0, np.random.RandomState(0)]
-    if NumpyVersion(np.__version__) >= '1.18.0':
-        rngs.append(np.random.default_rng(0))  # type: ignore
-    sizes = [(None, tuple()), (8, (8,)), ((4, 5, 6), (4, 5, 6))]
-
-    @pytest.mark.parametrize('rng', rngs)
-    @pytest.mark.parametrize('size_in, size_out', sizes)
-    def test_RVS(self, rng, size_in, size_out):
-        dist = stats.norm()
-        fni = stats.NumericalInverseHermite(dist)
-
-        rng2 = deepcopy(rng)
-        rvs = fni.rvs(size=size_in, random_state=rng)
-        assert(rvs.shape == size_out)
-
-        if rng2 is not None:
-            rng2 = check_random_state(rng2)
-            uniform = rng2.uniform(size=size_in)
-            rvs2 = stats.norm.ppf(uniform)
-            assert_allclose(rvs, rvs2)
-
-    if NumpyVersion(np.__version__) >= '1.18.0':
-        qrngs = [None, stats.qmc.Sobol(1, seed=0), stats.qmc.Halton(3, seed=0)]
-    else:
-        qrngs = []
-    # `size=None` should not add anything to the shape, `size=1` should
-    sizes = [(None, tuple()), (1, (1,)), (4, (4,)),
-             ((4,), (4,)), ((2, 4), (2, 4))]  # type: ignore
-    # Neither `d=None` nor `d=1` should add anything to the shape
-    ds = [(None, tuple()), (1, tuple()), (3, (3,))]
-
-    @pytest.mark.parametrize('qrng', qrngs)
-    @pytest.mark.parametrize('size_in, size_out', sizes)
-    @pytest.mark.parametrize('d_in, d_out', ds)
-    def test_QRVS(self, qrng, size_in, size_out, d_in, d_out):
-        dist = stats.norm()
-        fni = stats.NumericalInverseHermite(dist)
-
-        # If d and qrng.d are inconsistent, an error is raised
-        if d_in is not None and qrng is not None and qrng.d != d_in:
-            match = "`d` must be consistent with dimension of `qmc_engine`."
-            with pytest.raises(ValueError, match=match):
-                fni.qrvs(size_in, d=d_in, qmc_engine=qrng)
-            return
-
-        # Sometimes d is really determined by qrng
-        if d_in is None and qrng is not None and qrng.d != 1:
-            d_out = (qrng.d,)
-
-        shape_expected = size_out + d_out
-
-        qrng2 = deepcopy(qrng)
-        qrvs = fni.qrvs(size=size_in, d=d_in, qmc_engine=qrng)
-        assert(qrvs.shape == shape_expected)
-
-        if qrng2 is not None:
-            uniform = qrng2.random(np.prod(size_in) or 1)
-            qrvs2 = stats.norm.ppf(uniform).reshape(shape_expected)
-            assert_allclose(qrvs, qrvs2, atol=1e-12)
-
-    def test_QRVS_size_tuple(self):
-        # QMCEngine samples are always of shape (n, d). When `size` is a tuple,
-        # we set `n = prod(size)` in the call to qmc_engine.random, transform
-        # the sample, and reshape it to the final dimensions. When we reshape,
-        # we need to be careful, because the _columns_ of the sample returned
-        # by a QMCEngine are "independent"-ish, but the elements within the
-        # columns are not. We need to make sure that this doesn't get mixed up
-        # by reshaping: qrvs[..., i] should remain "independent"-ish of
-        # qrvs[..., i+1], but the elements within qrvs[..., i] should be
-        # transformed from the same low-discrepancy sequence.
-        if NumpyVersion(np.__version__) <= '1.18.0':
-            pytest.skip("QMC doesn't play well with old NumPy")
-
-        dist = stats.norm()
-        fni = stats.NumericalInverseHermite(dist)
-
-        size = (3, 4)
-        d = 5
-        qrng = stats.qmc.Halton(d, seed=0)
-        qrng2 = stats.qmc.Halton(d, seed=0)
-
-        uniform = qrng2.random(np.prod(size))
-
-        qrvs = fni.qrvs(size=size, d=d, qmc_engine=qrng)
-        qrvs2 = stats.norm.ppf(uniform)
-
-        for i in range(d):
-            sample = qrvs[..., i]
-            sample2 = qrvs2[:, i].reshape(size)
-            assert_allclose(sample, sample2, atol=1e-12)
-
-    def test_inaccurate_CDF(self):
-        # CDF function with inaccurate tail cannot be inverted; see gh-13319
-        # https://github.com/scipy/scipy/pull/13319#discussion_r626188955
-        shapes = (2.3098496451481823, 0.6268795430096368)
-        match = "The interpolating spline could not be created."
-
-        # fails with default tol
-        with pytest.raises(ValueError, match=match):
-            stats.NumericalInverseHermite(stats.beta(*shapes))
-
-        # no error with coarser tol
-        stats.NumericalInverseHermite(stats.beta(*shapes), tol=1e-10)
-
-    def test_custom_distribution(self):
-        class MyNormal:
-
-            def pdf(self, x):
-                return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
-
-            def cdf(self, x):
-                return special.ndtr(x)
-
-            def ppf(self, x):
-                return special.ndtri(x)
-
-        dist1 = MyNormal()
-        fni1 = stats.NumericalInverseHermite(dist1)
-
-        dist2 = stats.norm()
-        fni2 = stats.NumericalInverseHermite(dist2)
-
-        assert_allclose(fni1.rvs(random_state=0), fni2.rvs(random_state=0))
-
-
 class TestPageTrendTest:
     # expected statistic and p-values generated using R at
     # https://rdrr.io/cran/cultevo/, e.g.


Code review comments:
diff_hunk: @@ -0,0 +1,191 @@
+.. _sampling-hinv:
+
+Hermite interpolation based INVersion of CDF (HINV)
+===================================================
+
+.. currentmodule:: scipy.stats
+
+* Required: CDF
+* Optional: PDF, dPDF
+* Speed:
+
+  * Set-up: (very) slow
+  * Sampling: (very) fast
+
+HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+in each interval the inverse CDF is approximated by polynomials constructed by means
+of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+the accuracy by splitting a particular interval without recomputations in unaffected
+intervals. Three types of splines are implemented: linear, cubic, and quintic
+interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+also requires PDF and quintic interpolation PDF and its derivative.
+
+These splines have to be computed in a setup step. However, it only works for
+distributions with bounded domain; for distributions with unbounded domain the tails
+are chopped off such that the probability for the tail regions is small compared to
+the given u-resolution.
+
+The method is not exact, as it only produces random variates of the approximated
+distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
+quantile ``U`` i.e. ``X = approx_ppf(U)``) can be set to the
+required resolution (within machine precision). Notice that very small values of
+the u-resolution are possible but may increase the cost for the setup step.
+
+`NumericalInverseHermite` approximates the inverse of a continuous
+statistical distribution's CDF with a Hermite spline. Order of the
+hermite spline can be specified by passing the ``order`` parameter.
+
+As described in [1]_, it begins by evaluating the distribution's PDF and
+CDF at a mesh of quantiles ``x`` within the distribution's support.
+It uses the results to fit a Hermite spline ``H`` such that
+``H(p) == x``, where ``p`` is the array of percentiles corresponding
+with the quantiles ``x``. Therefore, the spline approximates the inverse
+of the distribution's CDF to machine precision at the percentiles ``p``,
+but typically, the spline will not be as accurate at the midpoints between
+the percentile points::
+
+    p_mid = (p[:-1] + p[1:])/2
+
+so the mesh of quantiles is refined as needed to reduce the maximum
+"u-error"::
+
+    u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
+
+below the specified tolerance `u_resolution`. Refinement stops when the required
+tolerance is achieved or when the number of mesh intervals after the next
+refinement could exceed the maximum allowed number `max_intervals`.
+
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+
+To create a generator to sample from the standard normal distribution, do:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+
+The `NumericalInverseHermite` has a method that approximates the PPF of the
+distribution.
+
+    >>> rng = NumericalInverseHermite(dist)
+    >>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
+    >>> np.allclose(rng.ppf(p), norm.ppf(p))
+    True
+
+Depending on the implementation of the distribution's random sampling
+method, the random variates generated may be nearly identical, given
+the same random state.
+
+    >>> dist = genexpon(9, 16, 3)
+    >>> rng = NumericalInverseHermite(dist)
+    >>> # `seed` ensures identical random streams are used by each `rvs` method
+    >>> seed = 500072020
+    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> np.allclose(rvs1, rvs2)
+    True
+
+To check that the random variates closely follow the given distribution, we can
+look at its histogram:
+
+.. plot::
+
+    >>> import numpy as np
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+    >>> import matplotlib.pyplot as plt
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+    >>> rvs = rng.rvs(10000)
+    >>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
+    >>> fx = norm.pdf(x)
+    >>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
+    >>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
+    >>> plt.xlabel('x')
+    >>> plt.ylabel('PDF(x)')
+    >>> plt.title('Numerical Inverse Hermite Samples')
+    >>> plt.legend()
+    >>> plt.show()
+
+Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use
+quintic Hermite interpolation to approximate the PPF by passing the ``order``
+parameter:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def dpdf(self, x):
+    ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
+
+Higher orders result is fewer number of intervals:
+
+    >>> rng3 = NumericalInverseHermite(dist, order=3)
+    >>> rng5 = NumericalInverseHermite(dist, order=5)
+    >>> rng3.intervals, rng5.intervals
+    (3000, 522)
+
+The u-error can be estimated by calling the `u_error` method. It runs a small
+Monte-Carlo simulation to estimate the u-error. By default, 100000 samples are
+used. To change this, you can pass the number of samples as an argument:
+
+    >>> rng1 = NumericalInverseHermite(dist, u_resolution=1.e-10)
+    >>> rng1.u_error(sample_size=1000000)  # uses one million samples
+    UError(max_error=9.53167544892608e-11, mean_absolute_error=2.2450136432146864e-11)
+
+This returns a namedtuple which contains the maximum u-error and the mean
+absolute u-error.
+
+The u-error can be reduced by decreasing the u-resolution (maximum allowed u-error):
+
+    >>> rng2 = NumericalInverseHermite(dist, u_resolution=1.e-13)
+    >>> rng2.u_error(sample_size=1000000)
+    UError(max_error=9.32027892364129e-14, mean_absolute_error=1.5194172675685075e-14)
+
+Note that this comes at the cost of increased setup time and number of intervals.
path: doc/source/tutorial/stats/sampling_hinv.rst
text: @user1:
Is this entirely correct? This phrasing implies that reducing the u-error by decreasing the u-resolution comes at the cost of the number of intervals, but they in fact increase, which seems to be a positive factor in reducing the u-error. 

I think you may have meant that 
```suggestion
Note that this comes at the cost of computation time as a result of the increased setup time and number of intervals. 
```

@author:
> Is this entirely correct? This phrasing implies that reducing the u-error by decreasing the u-resolution comes at the cost of the number of intervals, but they in fact increase, which seems to be a positive factor in reducing the u-error.

Ah, yes. This makes sense.

> I think you may have meant that

Yes, that's exactly what I meant!

@author:
Changed it this way, thanks!

diff_hunk: @@ -0,0 +1,191 @@
+.. _sampling-hinv:
+
+Hermite interpolation based INVersion of CDF (HINV)
+===================================================
+
+.. currentmodule:: scipy.stats
+
+* Required: CDF
+* Optional: PDF, dPDF
+* Speed:
+
+  * Set-up: (very) slow
+  * Sampling: (very) fast
+
+HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+in each interval the inverse CDF is approximated by polynomials constructed by means
+of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+the accuracy by splitting a particular interval without recomputations in unaffected
+intervals. Three types of splines are implemented: linear, cubic, and quintic
+interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+also requires PDF and quintic interpolation PDF and its derivative.
+
+These splines have to be computed in a setup step. However, it only works for
+distributions with bounded domain; for distributions with unbounded domain the tails
+are chopped off such that the probability for the tail regions is small compared to
+the given u-resolution.
+
+The method is not exact, as it only produces random variates of the approximated
+distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
+quantile ``U`` i.e. ``X = approx_ppf(U)``) can be set to the
+required resolution (within machine precision). Notice that very small values of
+the u-resolution are possible but may increase the cost for the setup step.
+
+`NumericalInverseHermite` approximates the inverse of a continuous
+statistical distribution's CDF with a Hermite spline. Order of the
+hermite spline can be specified by passing the ``order`` parameter.
+
+As described in [1]_, it begins by evaluating the distribution's PDF and
+CDF at a mesh of quantiles ``x`` within the distribution's support.
+It uses the results to fit a Hermite spline ``H`` such that
+``H(p) == x``, where ``p`` is the array of percentiles corresponding
+with the quantiles ``x``. Therefore, the spline approximates the inverse
+of the distribution's CDF to machine precision at the percentiles ``p``,
+but typically, the spline will not be as accurate at the midpoints between
+the percentile points::
+
+    p_mid = (p[:-1] + p[1:])/2
+
+so the mesh of quantiles is refined as needed to reduce the maximum
+"u-error"::
+
+    u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
+
+below the specified tolerance `u_resolution`. Refinement stops when the required
+tolerance is achieved or when the number of mesh intervals after the next
+refinement could exceed the maximum allowed number `max_intervals`.
+
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+
+To create a generator to sample from the standard normal distribution, do:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+
+The `NumericalInverseHermite` has a method that approximates the PPF of the
+distribution.
+
+    >>> rng = NumericalInverseHermite(dist)
+    >>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
+    >>> np.allclose(rng.ppf(p), norm.ppf(p))
+    True
+
+Depending on the implementation of the distribution's random sampling
+method, the random variates generated may be nearly identical, given
+the same random state.
+
+    >>> dist = genexpon(9, 16, 3)
+    >>> rng = NumericalInverseHermite(dist)
+    >>> # `seed` ensures identical random streams are used by each `rvs` method
+    >>> seed = 500072020
+    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> np.allclose(rvs1, rvs2)
+    True
+
+To check that the random variates closely follow the given distribution, we can
+look at its histogram:
+
+.. plot::
+
+    >>> import numpy as np
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+    >>> import matplotlib.pyplot as plt
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+    >>> rvs = rng.rvs(10000)
+    >>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
+    >>> fx = norm.pdf(x)
+    >>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
+    >>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
+    >>> plt.xlabel('x')
+    >>> plt.ylabel('PDF(x)')
+    >>> plt.title('Numerical Inverse Hermite Samples')
+    >>> plt.legend()
+    >>> plt.show()
+
+Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use
+quintic Hermite interpolation to approximate the PPF by passing the ``order``
+parameter:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def dpdf(self, x):
+    ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
+
+Higher orders result is fewer number of intervals:
+
+    >>> rng3 = NumericalInverseHermite(dist, order=3)
+    >>> rng5 = NumericalInverseHermite(dist, order=5)
+    >>> rng3.intervals, rng5.intervals
+    (3000, 522)
+
+The u-error can be estimated by calling the `u_error` method. It runs a small
+Monte-Carlo simulation to estimate the u-error. By default, 100000 samples are
+used. To change this, you can pass the number of samples as an argument:
+
+    >>> rng1 = NumericalInverseHermite(dist, u_resolution=1.e-10)
path: doc/source/tutorial/stats/sampling_hinv.rst
text: @user1:
Should the additional dot in `1.e-10` be consistent with `1e-10` like below? 

@author:
Yes, changed.

diff_hunk: @@ -1370,6 +1406,538 @@ cdef class NumericalInversePolynomial(Method):
         return UError(max_error, mae)
 
 
+cdef class NumericalInverseHermite(Method):
+    """
+    NumericalInverseHermite(dist, *, domain=None, order=3, u_resolution=1e-12, construction_points=None, max_intervals=100000, guide_factor=1, random_state=None)
+
+    Hermite interpolation based INVersion of CDF (HINV).
+
+    HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+    Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+    in each interval the inverse CDF is approximated by polynomials constructed by means
+    of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+    the accuracy by splitting a particular interval without recomputations in unaffected
+    intervals. Three types of splines are implemented: linear, cubic, and quintic
+    interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+    also requires PDF and quintic interpolation PDF and its derivative.
+
+    These splines have to be computed in a setup step. However, it only works for
+    distributions with bounded domain; for distributions with unbounded domain the tails
+    are chopped off such that the probability for the tail regions is small compared to
+    the given u-resolution.
+
+    The method is not exact, as it only produces random variates of the approximated
+    distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+    ``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
path: scipy/stats/_unuran/unuran_wrapper.pyx.templ
text: @user1:
```suggestion
    ``|U - CDF(X)|`` where ``X`` is the approximate percentile corresponding to the
```

@author:
Done.

diff_hunk: @@ -0,0 +1,191 @@
+.. _sampling-hinv:
+
+Hermite interpolation based INVersion of CDF (HINV)
+===================================================
+
+.. currentmodule:: scipy.stats
+
+* Required: CDF
+* Optional: PDF, dPDF
+* Speed:
+
+  * Set-up: (very) slow
+  * Sampling: (very) fast
+
+HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+in each interval the inverse CDF is approximated by polynomials constructed by means
+of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+the accuracy by splitting a particular interval without recomputations in unaffected
+intervals. Three types of splines are implemented: linear, cubic, and quintic
+interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+also requires PDF and quintic interpolation PDF and its derivative.
+
+These splines have to be computed in a setup step. However, it only works for
+distributions with bounded domain; for distributions with unbounded domain the tails
+are chopped off such that the probability for the tail regions is small compared to
+the given u-resolution.
+
+The method is not exact, as it only produces random variates of the approximated
+distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
path: doc/source/tutorial/stats/sampling_hinv.rst
text: @user1:
```suggestion
``|U - CDF(X)|`` where ``X`` is the approximate percentile corresponding to the
```

This is misspelled in the API doc as well

@author:
Corrected both here and in API docs.

diff_hunk: @@ -0,0 +1,191 @@
+.. _sampling-hinv:
+
+Hermite interpolation based INVersion of CDF (HINV)
+===================================================
+
+.. currentmodule:: scipy.stats
+
+* Required: CDF
+* Optional: PDF, dPDF
+* Speed:
+
+  * Set-up: (very) slow
+  * Sampling: (very) fast
+
+HINV is a variant of numerical inversion, where the inverse CDF is approximated using
+Hermite interpolation, i.e., the interval [0,1] is split into several intervals and
+in each interval the inverse CDF is approximated by polynomials constructed by means
+of values of the CDF and PDF at interval boundaries. This makes it possible to improve
+the accuracy by splitting a particular interval without recomputations in unaffected
+intervals. Three types of splines are implemented: linear, cubic, and quintic
+interpolation. For linear interpolation only the CDF is required. Cubic interpolation
+also requires PDF and quintic interpolation PDF and its derivative.
+
+These splines have to be computed in a setup step. However, it only works for
+distributions with bounded domain; for distributions with unbounded domain the tails
+are chopped off such that the probability for the tail regions is small compared to
+the given u-resolution.
+
+The method is not exact, as it only produces random variates of the approximated
+distribution. Nevertheless, the maximal numerical error in "u-direction" (i.e.
+``|U - CDF(X)|`` where ``X`` is the approximate percentile corressponding to the
+quantile ``U`` i.e. ``X = approx_ppf(U)``) can be set to the
+required resolution (within machine precision). Notice that very small values of
+the u-resolution are possible but may increase the cost for the setup step.
+
+`NumericalInverseHermite` approximates the inverse of a continuous
+statistical distribution's CDF with a Hermite spline. Order of the
+hermite spline can be specified by passing the ``order`` parameter.
+
+As described in [1]_, it begins by evaluating the distribution's PDF and
+CDF at a mesh of quantiles ``x`` within the distribution's support.
+It uses the results to fit a Hermite spline ``H`` such that
+``H(p) == x``, where ``p`` is the array of percentiles corresponding
+with the quantiles ``x``. Therefore, the spline approximates the inverse
+of the distribution's CDF to machine precision at the percentiles ``p``,
+but typically, the spline will not be as accurate at the midpoints between
+the percentile points::
+
+    p_mid = (p[:-1] + p[1:])/2
+
+so the mesh of quantiles is refined as needed to reduce the maximum
+"u-error"::
+
+    u_error = np.max(np.abs(dist.cdf(H(p_mid)) - p_mid))
+
+below the specified tolerance `u_resolution`. Refinement stops when the required
+tolerance is achieved or when the number of mesh intervals after the next
+refinement could exceed the maximum allowed number `max_intervals`.
+
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+
+To create a generator to sample from the standard normal distribution, do:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+
+The `NumericalInverseHermite` has a method that approximates the PPF of the
+distribution.
+
+    >>> rng = NumericalInverseHermite(dist)
+    >>> p = np.linspace(0.01, 0.99, 99) # percentiles from 1% to 99%
+    >>> np.allclose(rng.ppf(p), norm.ppf(p))
+    True
+
+Depending on the implementation of the distribution's random sampling
+method, the random variates generated may be nearly identical, given
+the same random state.
+
+    >>> dist = genexpon(9, 16, 3)
+    >>> rng = NumericalInverseHermite(dist)
+    >>> # `seed` ensures identical random streams are used by each `rvs` method
+    >>> seed = 500072020
+    >>> rvs1 = dist.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> rvs2 = rng.rvs(size=100, random_state=np.random.default_rng(seed))
+    >>> np.allclose(rvs1, rvs2)
+    True
+
+To check that the random variates closely follow the given distribution, we can
+look at its histogram:
+
+.. plot::
+
+    >>> import numpy as np
+    >>> from scipy.stats import NumericalInverseHermite
+    >>> from scipy.stats import norm, genexpon
+    >>> from scipy.special import ndtr
+    >>> import matplotlib.pyplot as plt
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, random_state=urng)
+    >>> rvs = rng.rvs(10000)
+    >>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, 1000)
+    >>> fx = norm.pdf(x)
+    >>> plt.plot(x, fx, 'r-', lw=2, label='true distribution')
+    >>> plt.hist(rvs, bins=20, density=True, alpha=0.8, label='random variates')
+    >>> plt.xlabel('x')
+    >>> plt.ylabel('PDF(x)')
+    >>> plt.title('Numerical Inverse Hermite Samples')
+    >>> plt.legend()
+    >>> plt.show()
+
+Given the derivative of the PDF w.r.t the variate (i.e. ``x``), we can use
+quintic Hermite interpolation to approximate the PPF by passing the ``order``
+parameter:
+
+    >>> class StandardNormal:
+    ...     def pdf(self, x):
+    ...        return 1/np.sqrt(2*np.pi) * np.exp(-x**2 / 2)
+    ...     def dpdf(self, x):
+    ...        return -1/np.sqrt(2*np.pi) * x * np.exp(-x**2 / 2)
+    ...     def cdf(self, x):
+    ...        return ndtr(x)
+    ... 
+    >>> dist = StandardNormal()
+    >>> urng = np.random.default_rng()
+    >>> rng = NumericalInverseHermite(dist, order=5, random_state=urng)
+
+Higher orders result is fewer number of intervals:
path: doc/source/tutorial/stats/sampling_hinv.rst
text: @user1:
```suggestion
Higher orders result in a fewer number of intervals:
```

very nitpicky suggestion but "result" is an intransitive verb in this usage 

@author:
Done.