image-match/python/py/Lib/site-packages/scipy/stats/tests/test_sampling.py

import threading
import pickle
import pytest
from copy import deepcopy
import platform
import sys
import math
import warnings

import numpy as np
from numpy.testing import assert_allclose, assert_equal
from scipy.stats.sampling import (
    TransformedDensityRejection,
    DiscreteAliasUrn,
    DiscreteGuideTable,
    NumericalInversePolynomial,
    NumericalInverseHermite,
    RatioUniforms,
    SimpleRatioUniforms,
    UNURANError
)
from pytest import raises as assert_raises
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.


# Normal distribution shared between all the continuous methods
class StandardNormal:
    def pdf(self, x):
        # normalization constant needed for NumericalInverseHermite
        return 1./np.sqrt(2.*np.pi) * np.exp(-0.5 * x*x)

    def dpdf(self, x):
        return 1./np.sqrt(2.*np.pi) * -x * np.exp(-0.5 * x*x)

    def cdf(self, x):
        return special.ndtr(x)


all_methods = [
    ("TransformedDensityRejection", {"dist": StandardNormal()}),
    ("DiscreteAliasUrn", {"dist": [0.02, 0.18, 0.8]}),
    ("DiscreteGuideTable", {"dist": [0.02, 0.18, 0.8]}),
    ("NumericalInversePolynomial", {"dist": StandardNormal()}),
    ("NumericalInverseHermite", {"dist": StandardNormal()}),
    ("SimpleRatioUniforms", {"dist": StandardNormal(), "mode": 0})
]

if (sys.implementation.name == 'pypy'
        and sys.implementation.version < (7, 3, 10)):
    # changed in PyPy for v7.3.10
    floaterr = r"unsupported operand type for float\(\): 'list'"
else:
    floaterr = r"must be real number, not list"
# Make sure an internal error occurs in UNU.RAN when invalid callbacks are
# passed. Moreover, different generators throw different error messages.
# So, in case of an `UNURANError`, we do not validate the error message.
bad_pdfs_common = [
    # Negative PDF
    (lambda x: -x, UNURANError, r"..."),
    # Returning wrong type
    (lambda x: [], TypeError, floaterr),
    # Undefined name inside the function
    (lambda x: foo, NameError, r"name 'foo' is not defined"),  # type: ignore[name-defined]  # noqa: F821, E501
    # Infinite value returned => Overflow error.
    (lambda x: np.inf, UNURANError, r"..."),
    # NaN value => internal error in UNU.RAN
    (lambda x: np.nan, UNURANError, r"..."),
    # signature of PDF wrong
    (lambda: 1.0, TypeError, r"takes 0 positional arguments but 1 was given")
]


# same approach for dpdf
bad_dpdf_common = [
    # Infinite value returned.
    (lambda x: np.inf, UNURANError, r"..."),
    # NaN value => internal error in UNU.RAN
    (lambda x: np.nan, UNURANError, r"..."),
    # Returning wrong type
    (lambda x: [], TypeError, floaterr),
    # Undefined name inside the function
    (lambda x: foo, NameError, r"name 'foo' is not defined"),  # type: ignore[name-defined]  # noqa: F821, E501
    # signature of dPDF wrong
    (lambda: 1.0, TypeError, r"takes 0 positional arguments but 1 was given")
]


# same approach for logpdf
bad_logpdfs_common = [
    # Returning wrong type
    (lambda x: [], TypeError, floaterr),
    # Undefined name inside the function
    (lambda x: foo, NameError, r"name 'foo' is not defined"),  # type: ignore[name-defined]  # noqa: F821, E501
    # Infinite value returned => Overflow error.
    (lambda x: np.inf, UNURANError, r"..."),
    # NaN value => internal error in UNU.RAN
    (lambda x: np.nan, UNURANError, r"..."),
    # signature of logpdf wrong
    (lambda: 1.0, TypeError, r"takes 0 positional arguments but 1 was given")
]


bad_pv_common = [
    ([], r"must contain at least one element"),
    ([[1.0, 0.0]], r"wrong number of dimensions \(expected 1, got 2\)"),
    ([0.2, 0.4, np.nan, 0.8], r"must contain only finite / non-nan values"),
    ([0.2, 0.4, np.inf, 0.8], r"must contain only finite / non-nan values"),
    ([0.0, 0.0], r"must contain at least one non-zero value"),
]


# size of the domains is incorrect
bad_sized_domains = [
    # > 2 elements in the domain
    ((1, 2, 3), ValueError, r"must be a length 2 tuple"),
    # empty domain
    ((), ValueError, r"must be a length 2 tuple")
]

# domain values are incorrect
bad_domains = [
    ((2, 1), UNURANError, r"left >= right"),
    ((1, 1), UNURANError, r"left >= right"),
]

# infinite and nan values present in domain.
inf_nan_domains = [
    # left >= right
    ((10, 10), UNURANError, r"left >= right"),
    ((np.inf, np.inf), UNURANError, r"left >= right"),
    ((-np.inf, -np.inf), UNURANError, r"left >= right"),
    ((np.inf, -np.inf), UNURANError, r"left >= right"),
    # Also include nans in some of the domains.
    ((-np.inf, np.nan), ValueError, r"only non-nan values"),
    ((np.nan, np.inf), ValueError, r"only non-nan values")
]

# `nan` values present in domain. Some distributions don't support
# infinite tails, so don't mix the nan values with infinities.
nan_domains = [
    ((0, np.nan), ValueError, r"only non-nan values"),
    ((np.nan, np.nan), ValueError, r"only non-nan values")
]


# all the methods should throw errors for nan, bad sized, and bad valued
# domains.
@pytest.mark.parametrize("domain, err, msg",
                         bad_domains + bad_sized_domains +
                         nan_domains)  # type: ignore[operator]
@pytest.mark.parametrize("method, kwargs", all_methods)
def test_bad_domain(domain, err, msg, method, kwargs):
    Method = getattr(stats.sampling, method)
    with pytest.raises(err, match=msg):
        Method(**kwargs, domain=domain)


@pytest.mark.parametrize("method, kwargs", all_methods)
def test_random_state(method, kwargs):
    Method = getattr(stats.sampling, method)

    # simple seed that works for any version of NumPy
    seed = 123
    rng1 = Method(**kwargs, random_state=seed)
    rng2 = Method(**kwargs, random_state=seed)
    assert_equal(rng1.rvs(100), rng2.rvs(100))

    # global seed
    rng = np.random.RandomState(123)
    rng1 = Method(**kwargs)
    rvs1 = rng1.rvs(100, random_state=rng)
    np.random.seed(None)  # valid use of np.random.seed
    rng2 = Method(**kwargs, random_state=123)
    rvs2 = rng2.rvs(100)
    assert_equal(rvs1, rvs2)

    # Generator seed for new NumPy
    # when a RandomState is given, it should take the bitgen_t
    # member of the class and create a Generator instance.
    seed1 = np.random.RandomState(np.random.MT19937(123))
    seed2 = np.random.Generator(np.random.MT19937(123))
    rng1 = Method(**kwargs, random_state=seed1)
    rng2 = Method(**kwargs, random_state=seed2)
    assert_equal(rng1.rvs(100), rng2.rvs(100))


def test_set_random_state():
    rng1 = TransformedDensityRejection(StandardNormal(), random_state=123)
    rng2 = TransformedDensityRejection(StandardNormal())
    rng2.set_random_state(123)
    assert_equal(rng1.rvs(100), rng2.rvs(100))
    rng = TransformedDensityRejection(StandardNormal(), random_state=123)
    rvs1 = rng.rvs(100)
    rng.set_random_state(123)
    rvs2 = rng.rvs(100)
    assert_equal(rvs1, rvs2)


def test_threading_behaviour():
    # Test if the API is thread-safe.
    # This verifies if the lock mechanism and the use of `PyErr_Occurred`
    # is correct.
    errors = {"err1": None, "err2": None}

    class Distribution:
        def __init__(self, pdf_msg):
            self.pdf_msg = pdf_msg

        def pdf(self, x):
            if 49.9 < x < 50.0:
                raise ValueError(self.pdf_msg)
            return x

        def dpdf(self, x):
            return 1

    def func1():
        dist = Distribution('foo')
        rng = TransformedDensityRejection(dist, domain=(10, 100),
                                          random_state=12)
        try:
            rng.rvs(100000)
        except ValueError as e:
            errors['err1'] = e.args[0]

    def func2():
        dist = Distribution('bar')
        rng = TransformedDensityRejection(dist, domain=(10, 100),
                                          random_state=2)
        try:
            rng.rvs(100000)
        except ValueError as e:
            errors['err2'] = e.args[0]

    t1 = threading.Thread(target=func1)
    t2 = threading.Thread(target=func2)

    t1.start()
    t2.start()

    t1.join()
    t2.join()

    assert errors['err1'] == 'foo'
    assert errors['err2'] == 'bar'


@pytest.mark.parametrize("method, kwargs", all_methods)
def test_pickle(method, kwargs):
    Method = getattr(stats.sampling, method)
    rng1 = Method(**kwargs, random_state=123)
    obj = pickle.dumps(rng1)
    rng2 = pickle.loads(obj)
    assert_equal(rng1.rvs(100), rng2.rvs(100))


@pytest.mark.parametrize("size", [None, 0, (0, ), 1, (10, 3), (2, 3, 4, 5),
                                  (0, 0), (0, 1)])
def test_rvs_size(size):
    # As the `rvs` method is present in the base class and shared between
    # all the classes, we can just test with one of the methods.
    rng = TransformedDensityRejection(StandardNormal())
    if size is None:
        assert np.isscalar(rng.rvs(size))
    else:
        if np.isscalar(size):
            size = (size, )
        assert rng.rvs(size).shape == size


def test_with_scipy_distribution():
    # test if the setup works with SciPy's rv_frozen distributions
    dist = stats.norm()
    urng = np.random.default_rng(0)
    rng = NumericalInverseHermite(dist, random_state=urng)
    u = np.linspace(0, 1, num=100)
    check_cont_samples(rng, dist, dist.stats())
    assert_allclose(dist.ppf(u), rng.ppf(u))
    # test if it works with `loc` and `scale`
    dist = stats.norm(loc=10., scale=5.)
    rng = NumericalInverseHermite(dist, random_state=urng)
    check_cont_samples(rng, dist, dist.stats())
    assert_allclose(dist.ppf(u), rng.ppf(u))
    # check for discrete distributions
    dist = stats.binom(10, 0.2)
    rng = DiscreteAliasUrn(dist, random_state=urng)
    domain = dist.support()
    pv = dist.pmf(np.arange(domain[0], domain[1]+1))
    check_discr_samples(rng, pv, dist.stats())


def check_cont_samples(rng, dist, mv_ex, rtol=1e-7, atol=1e-1):
    rvs = rng.rvs(100000)
    mv = rvs.mean(), rvs.var()
    # test the moments only if the variance is finite
    if np.isfinite(mv_ex[1]):
        assert_allclose(mv, mv_ex, rtol=rtol, atol=atol)
    # Cramer Von Mises test for goodness-of-fit
    rvs = rng.rvs(500)
    dist.cdf = np.vectorize(dist.cdf)
    pval = cramervonmises(rvs, dist.cdf).pvalue
    assert pval > 0.1


def check_discr_samples(rng, pv, mv_ex, rtol=1e-3, atol=1e-1):
    rvs = rng.rvs(100000)
    # test if the first few moments match
    mv = rvs.mean(), rvs.var()
    assert_allclose(mv, mv_ex, rtol=rtol, atol=atol)
    # normalize
    pv = pv / pv.sum()
    # chi-squared test for goodness-of-fit
    obs_freqs = np.zeros_like(pv)
    _, freqs = np.unique(rvs, return_counts=True)
    freqs = freqs / freqs.sum()
    obs_freqs[:freqs.size] = freqs
    pval = chisquare(obs_freqs, pv).pvalue
    assert pval > 0.1


def test_warning_center_not_in_domain():
    # UNURAN will warn if the center provided or the one computed w/o the
    # domain is outside of the domain
    msg = "102 : center moved into domain of distribution"
    with pytest.warns(RuntimeWarning, match=msg):
        NumericalInversePolynomial(StandardNormal(), center=0, domain=(3, 5))
    with pytest.warns(RuntimeWarning, match=msg):
        NumericalInversePolynomial(StandardNormal(), domain=(3, 5))


@pytest.mark.parametrize('method', ["SimpleRatioUniforms",
                                    "NumericalInversePolynomial",
                                    "TransformedDensityRejection"])
def test_error_mode_not_in_domain(method):
    # UNURAN raises an error if the mode is not in the domain
    # the behavior is different compared to the case that center is not in the
    # domain. mode is supposed to be the exact value, center can be an
    # approximate value
    Method = getattr(stats.sampling, method)
    msg = "17 : mode not in domain"
    with pytest.raises(UNURANError, match=msg):
        Method(StandardNormal(), mode=0, domain=(3, 5))


@pytest.mark.parametrize('method', ["NumericalInverseHermite",
                                    "NumericalInversePolynomial"])
class TestQRVS:
    def test_input_validation(self, method):
        match = "`qmc_engine` must be an instance of..."
        with pytest.raises(ValueError, match=match):
            Method = getattr(stats.sampling, method)
            gen = Method(StandardNormal())
            gen.qrvs(qmc_engine=0)

        # issues with QMCEngines and old NumPy
        Method = getattr(stats.sampling, method)
        gen = Method(StandardNormal())

        match = "`d` must be consistent with dimension of `qmc_engine`."
        with pytest.raises(ValueError, match=match):
            gen.qrvs(d=3, qmc_engine=stats.qmc.Halton(2))

    qrngs = [None, stats.qmc.Sobol(1, seed=0), stats.qmc.Halton(3, seed=0)]
    # `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)
    @pytest.mark.thread_unsafe(reason="fails in parallel")
    def test_QRVS_shape_consistency(self, qrng, size_in, size_out,
                                    d_in, d_out, method):
        w32 = sys.platform == "win32" and platform.architecture()[0] == "32bit"
        if w32 and method == "NumericalInversePolynomial":
            pytest.xfail("NumericalInversePolynomial.qrvs fails for Win "
                         "32-bit")

        dist = StandardNormal()
        Method = getattr(stats.sampling, method)
        gen = Method(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):
                gen.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 = gen.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, method):
        # 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.

        dist = StandardNormal()
        Method = getattr(stats.sampling, method)
        gen = Method(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 = gen.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)


class TestTransformedDensityRejection:
    # Simple Custom Distribution
    class dist0:
        def pdf(self, x):
            return 3/4 * (1-x*x)

        def dpdf(self, x):
            return 3/4 * (-2*x)

        def cdf(self, x):
            return 3/4 * (x - x**3/3 + 2/3)

        def support(self):
            return -1, 1

    # Standard Normal Distribution
    class dist1:
        def pdf(self, x):
            return stats.norm._pdf(x / 0.1)

        def dpdf(self, x):
            return -x / 0.01 * stats.norm._pdf(x / 0.1)

        def cdf(self, x):
            return stats.norm._cdf(x / 0.1)

    # pdf with piecewise linear function as transformed density
    # with T = -1/sqrt with shift. Taken from UNU.RAN test suite
    # (from file t_tdr_ps.c)
    class dist2:
        def __init__(self, shift):
            self.shift = shift

        def pdf(self, x):
            x -= self.shift
            y = 1. / (abs(x) + 1.)
            return 0.5 * y * y

        def dpdf(self, x):
            x -= self.shift
            y = 1. / (abs(x) + 1.)
            y = y * y * y
            return y if (x < 0.) else -y

        def cdf(self, x):
            x -= self.shift
            if x <= 0.:
                return 0.5 / (1. - x)
            else:
                return 1. - 0.5 / (1. + x)

    dists = [dist0(), dist1(), dist2(0.), dist2(10000.)]

    # exact mean and variance of the distributions in the list dists
    mv0 = [0., 4./15.]
    mv1 = [0., 0.01]
    mv2 = [0., np.inf]
    mv3 = [10000., np.inf]
    mvs = [mv0, mv1, mv2, mv3]

    @pytest.mark.parametrize("dist, mv_ex",
                             zip(dists, mvs))
    @pytest.mark.thread_unsafe(reason="deadlocks for unknown reasons")
    def test_basic(self, dist, mv_ex):
        with warnings.catch_warnings():
            # filter the warnings thrown by UNU.RAN
            warnings.simplefilter("ignore", RuntimeWarning)
            rng = TransformedDensityRejection(dist, random_state=42)
        check_cont_samples(rng, dist, mv_ex)

    # PDF 0 everywhere => bad construction points
    bad_pdfs = [(lambda x: 0, UNURANError, r"50 : bad construction points.")]
    bad_pdfs += bad_pdfs_common  # type: ignore[arg-type]

    @pytest.mark.parametrize("pdf, err, msg", bad_pdfs)
    def test_bad_pdf(self, pdf, err, msg):
        class dist:
            pass
        dist.pdf = pdf
        dist.dpdf = lambda x: 1  # an arbitrary dPDF
        with pytest.raises(err, match=msg):
            TransformedDensityRejection(dist)

    @pytest.mark.parametrize("dpdf, err, msg", bad_dpdf_common)
    def test_bad_dpdf(self, dpdf, err, msg):
        class dist:
            pass
        dist.pdf = lambda x: x
        dist.dpdf = dpdf
        with pytest.raises(err, match=msg):
            TransformedDensityRejection(dist, domain=(1, 10))

    # 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):
            TransformedDensityRejection(StandardNormal(), domain=domain)

    @pytest.mark.parametrize("construction_points", [-1, 0, 0.1])
    def test_bad_construction_points_scalar(self, construction_points):
        with pytest.raises(ValueError, match=r"`construction_points` must be "
                                             r"a positive integer."):
            TransformedDensityRejection(
                StandardNormal(), construction_points=construction_points
            )

    def test_bad_construction_points_array(self):
        # empty array
        construction_points = []
        with pytest.raises(ValueError, match=r"`construction_points` must "
                                             r"either be a "
                                             r"scalar or a non-empty array."):
            TransformedDensityRejection(
                StandardNormal(), construction_points=construction_points
            )

        # construction_points not monotonically increasing
        construction_points = [1, 1, 1, 1, 1, 1]
        with pytest.warns(RuntimeWarning, match=r"33 : starting points not "
                                                r"strictly monotonically "
                                                r"increasing"):
            TransformedDensityRejection(
                StandardNormal(), construction_points=construction_points
            )

        # construction_points containing nans
        construction_points = [np.nan, np.nan, np.nan]
        with pytest.raises(UNURANError, match=r"50 : bad construction "
                                              r"points."):
            TransformedDensityRejection(
                StandardNormal(), construction_points=construction_points
            )

        # construction_points out of domain
        construction_points = [-10, 10]
        with pytest.warns(RuntimeWarning, match=r"50 : starting point out of "
                                                r"domain"):
            TransformedDensityRejection(
                StandardNormal(), domain=(-3, 3),
                construction_points=construction_points
            )

    @pytest.mark.parametrize("c", [-1., np.nan, np.inf, 0.1, 1.])
    def test_bad_c(self, c):
        msg = r"`c` must either be -0.5 or 0."
        with pytest.raises(ValueError, match=msg):
            TransformedDensityRejection(StandardNormal(), c=-1.)

    u = [np.linspace(0, 1, num=1000), [], [[]], [np.nan],
         [-np.inf, np.nan, np.inf], 0,
         [[np.nan, 0.5, 0.1], [0.2, 0.4, np.inf], [-2, 3, 4]]]

    @pytest.mark.parametrize("u", u)
    def test_ppf_hat(self, u):
        # Increase the `max_squeeze_hat_ratio` so the ppf_hat is more
        # accurate.
        rng = TransformedDensityRejection(StandardNormal(),
                                          max_squeeze_hat_ratio=0.9999)
        # Older versions of NumPy throw RuntimeWarnings for comparisons
        # with nan.
        with warnings.catch_warnings():
            msg = "invalid value encountered in "
            warnings.filterwarnings("ignore", msg + "greater", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "greater_equal", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less_equal", RuntimeWarning)
            res = rng.ppf_hat(u)
            expected = stats.norm.ppf(u)
        assert_allclose(res, expected, rtol=1e-3, atol=1e-5)
        assert res.shape == expected.shape

    def test_bad_dist(self):
        # Empty distribution
        class dist:
            ...

        msg = r"`pdf` required but not found."
        with pytest.raises(ValueError, match=msg):
            TransformedDensityRejection(dist)

        # dPDF not present in dist
        class dist:
            pdf = lambda x: 1-x*x  # noqa: E731

        msg = r"`dpdf` required but not found."
        with pytest.raises(ValueError, match=msg):
            TransformedDensityRejection(dist)


class TestDiscreteAliasUrn:
    # DAU fails on these probably because of large domains and small
    # computation errors in PMF. Mean/SD match but chi-squared test fails.
    basic_fail_dists = {
        'nchypergeom_fisher',  # numerical errors on tails
        'nchypergeom_wallenius',  # numerical errors on tails
        'randint'  # fails on 32-bit ubuntu
    }

    @pytest.mark.parametrize("distname, params", distdiscrete)
    def test_basic(self, distname, params):
        if distname in self.basic_fail_dists:
            msg = ("DAU fails on these probably because of large domains "
                   "and small computation errors in PMF.")
            pytest.skip(msg)
        if not isinstance(distname, str):
            dist = distname
        else:
            dist = getattr(stats, distname)
        dist = dist(*params)
        domain = dist.support()
        if not np.isfinite(domain[1] - domain[0]):
            # DAU only works with finite domain. So, skip the distributions
            # with infinite tails.
            pytest.skip("DAU only works with a finite domain.")
        k = np.arange(domain[0], domain[1]+1)
        pv = dist.pmf(k)
        mv_ex = dist.stats('mv')
        rng = DiscreteAliasUrn(dist, random_state=42)
        check_discr_samples(rng, pv, mv_ex)

    # Can't use bad_pmf_common here as we evaluate PMF early on to avoid
    # unhelpful errors from UNU.RAN.
    bad_pmf = [
        # inf returned
        (lambda x: np.inf, ValueError,
         r"must contain only finite / non-nan values"),
        # nan returned
        (lambda x: np.nan, ValueError,
         r"must contain only finite / non-nan values"),
        # all zeros
        (lambda x: 0.0, ValueError,
         r"must contain at least one non-zero value"),
        # Undefined name inside the function
        (lambda x: foo, NameError,  # type: ignore[name-defined]  # noqa: F821
         r"name 'foo' is not defined"),
        # Returning wrong type.
        (lambda x: [], ValueError,
         r"setting an array element with a sequence."),
        # probabilities < 0
        (lambda x: -x, UNURANError,
         r"50 : probability < 0"),
        # signature of PMF wrong
        (lambda: 1.0, TypeError,
         r"takes 0 positional arguments but 1 was given")
    ]

    @pytest.mark.parametrize("pmf, err, msg", bad_pmf)
    def test_bad_pmf(self, pmf, err, msg):
        class dist:
            pass
        dist.pmf = pmf
        with pytest.raises(err, match=msg):
            DiscreteAliasUrn(dist, domain=(1, 10))

    @pytest.mark.parametrize("pv", [[0.18, 0.02, 0.8],
                                    [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]])
    def test_sampling_with_pv(self, pv):
        pv = np.asarray(pv, dtype=np.float64)
        rng = DiscreteAliasUrn(pv, random_state=123)
        rng.rvs(100_000)
        pv = pv / pv.sum()
        variates = np.arange(0, len(pv))
        # test if the first few moments match
        m_expected = np.average(variates, weights=pv)
        v_expected = np.average((variates - m_expected) ** 2, weights=pv)
        mv_expected = m_expected, v_expected
        check_discr_samples(rng, pv, mv_expected)

    @pytest.mark.parametrize("pv, msg", bad_pv_common)
    def test_bad_pv(self, pv, msg):
        with pytest.raises(ValueError, match=msg):
            DiscreteAliasUrn(pv)

    # DAU doesn't support infinite tails. So, it should throw an error when
    # inf is present in the domain.
    inf_domain = [(-np.inf, np.inf), (np.inf, np.inf), (-np.inf, -np.inf),
                  (0, np.inf), (-np.inf, 0)]

    @pytest.mark.parametrize("domain", inf_domain)
    def test_inf_domain(self, domain):
        with pytest.raises(ValueError, match=r"must be finite"):
            DiscreteAliasUrn(stats.binom(10, 0.2), domain=domain)

    def test_bad_urn_factor(self):
        with pytest.warns(RuntimeWarning, match=r"relative urn size < 1."):
            DiscreteAliasUrn([0.5, 0.5], urn_factor=-1)

    def test_bad_args(self):
        msg = (r"`domain` must be provided when the "
               r"probability vector is not available.")

        class dist:
            def pmf(self, x):
                return x

        with pytest.raises(ValueError, match=msg):
            DiscreteAliasUrn(dist)

    def test_gh19359(self):
        pv = special.softmax(np.ones((1533,)))
        rng = DiscreteAliasUrn(pv, random_state=42)
        # check the correctness
        check_discr_samples(rng, pv, (1532 / 2, (1532**2 - 1) / 12),
                            rtol=5e-3)


class TestNumericalInversePolynomial:
    # Simple Custom Distribution
    class dist0:
        def pdf(self, x):
            return 3/4 * (1-x*x)

        def cdf(self, x):
            return 3/4 * (x - x**3/3 + 2/3)

        def support(self):
            return -1, 1

    # Standard Normal Distribution
    class dist1:
        def pdf(self, x):
            return stats.norm._pdf(x / 0.1)

        def cdf(self, x):
            return stats.norm._cdf(x / 0.1)

    # Sin 2 distribution
    #          /  0.05 + 0.45*(1 +sin(2 Pi x))  if |x| <= 1
    #  f(x) = <
    #          \  0        otherwise
    # Taken from UNU.RAN test suite (from file t_pinv.c)
    class dist2:
        def pdf(self, x):
            return 0.05 + 0.45 * (1 + np.sin(2*np.pi*x))

        def cdf(self, x):
            return (0.05*(x + 1) +
                    0.9*(1. + 2.*np.pi*(1 + x) - np.cos(2.*np.pi*x)) /
                    (4.*np.pi))

        def support(self):
            return -1, 1

    # Sin 10 distribution
    #          /  0.05 + 0.45*(1 +sin(2 Pi x))  if |x| <= 5
    #  f(x) = <
    #          \  0        otherwise
    # Taken from UNU.RAN test suite (from file t_pinv.c)
    class dist3:
        def pdf(self, x):
            return 0.2 * (0.05 + 0.45 * (1 + np.sin(2*np.pi*x)))

        def cdf(self, x):
            return x/10. + 0.5 + 0.09/(2*np.pi) * (np.cos(10*np.pi) -
                                                   np.cos(2*np.pi*x))

        def support(self):
            return -5, 5

    dists = [dist0(), dist1(), dist2(), dist3()]

    # exact mean and variance of the distributions in the list dists
    mv0 = [0., 4./15.]
    mv1 = [0., 0.01]
    mv2 = [-0.45/np.pi, 2/3*0.5 - 0.45**2/np.pi**2]
    mv3 = [-0.45/np.pi, 0.2 * 250/3 * 0.5 - 0.45**2/np.pi**2]
    mvs = [mv0, mv1, mv2, mv3]

    @pytest.mark.thread_unsafe(reason="deadlocks for unknown reasons")
    @pytest.mark.parametrize("dist, mv_ex",
                             zip(dists, mvs))
    def test_basic(self, dist, mv_ex):
        rng = NumericalInversePolynomial(dist, random_state=42)
        check_cont_samples(rng, dist, mv_ex)

    @pytest.mark.xslow
    @pytest.mark.parametrize("distname, params", distcont)
    def test_basic_all_scipy_dists(self, distname, params):

        very_slow_dists = ['anglit', 'gausshyper', 'kappa4',
                           'ksone', 'kstwo', 'levy_l',
                           'levy_stable', 'studentized_range',
                           'trapezoid', 'triang', 'vonmises']
        # for these distributions, some assertions fail due to minor
        # numerical differences. They can be avoided either by changing
        # the seed or by increasing the u_resolution.
        fail_dists = ['chi2', 'fatiguelife', 'gibrat',
                      'halfgennorm', 'lognorm', 'ncf',
                      'ncx2', 'pareto', 't']
        # for these distributions, skip the check for agreement between sample
        # moments and true moments. We cannot expect them to pass due to the
        # high variance of sample moments.
        skip_sample_moment_check = ['rel_breitwigner']

        if distname in very_slow_dists:
            pytest.skip(f"PINV too slow for {distname}")
        if distname in fail_dists:
            pytest.skip(f"PINV fails for {distname}")
        dist = (getattr(stats, distname)
                if isinstance(distname, str)
                else distname)
        dist = dist(*params)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)
            rng = NumericalInversePolynomial(dist, random_state=42)
        if distname in skip_sample_moment_check:
            return
        check_cont_samples(rng, dist, [dist.mean(), dist.var()])

    @pytest.mark.parametrize("pdf, err, msg", bad_pdfs_common)
    def test_bad_pdf(self, pdf, err, msg):
        class dist:
            pass
        dist.pdf = pdf
        with pytest.raises(err, match=msg):
            NumericalInversePolynomial(dist, domain=[0, 5])

    @pytest.mark.parametrize("logpdf, err, msg", bad_logpdfs_common)
    def test_bad_logpdf(self, logpdf, err, msg):
        class dist:
            pass
        dist.logpdf = logpdf
        with pytest.raises(err, match=msg):
            NumericalInversePolynomial(dist, domain=[0, 5])

    # 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):
            NumericalInversePolynomial(StandardNormal(), domain=domain)

    u = [
        # test if quantile 0 and 1 return -inf and inf respectively and check
        # the correctness of the PPF for equidistant points between 0 and 1.
        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 = NumericalInversePolynomial(dist, u_resolution=1e-14)
        # Older versions of NumPy throw RuntimeWarnings for comparisons
        # with nan.
        with warnings.catch_warnings():
            msg = "invalid value encountered in "
            warnings.filterwarnings("ignore", msg + "greater", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "greater_equal", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less_equal", RuntimeWarning)
            res = rng.ppf(u)
            expected = stats.norm.ppf(u)
        assert_allclose(res, expected, rtol=1e-11, atol=1e-11)
        assert res.shape == expected.shape

    x = [np.linspace(-10, 10, num=10000), [], [[]], [np.nan],
         [-np.inf, np.nan, np.inf], 0,
         [[np.nan, 0.5, 0.1], [0.2, 0.4, np.inf], [-np.inf, 3, 4]]]

    @pytest.mark.parametrize("x", x)
    def test_cdf(self, x):
        dist = StandardNormal()
        rng = NumericalInversePolynomial(dist, u_resolution=1e-14)
        # Older versions of NumPy throw RuntimeWarnings for comparisons
        # with nan.
        with warnings.catch_warnings():
            msg = "invalid value encountered in "
            warnings.filterwarnings("ignore", msg + "greater", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "greater_equal", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less_equal", RuntimeWarning)
            res = rng.cdf(x)
            expected = stats.norm.cdf(x)
        assert_allclose(res, expected, rtol=1e-11, atol=1e-11)
        assert res.shape == expected.shape

    @pytest.mark.slow
    def test_u_error(self):
        dist = StandardNormal()
        rng = NumericalInversePolynomial(dist, u_resolution=1e-10)
        max_error, mae = rng.u_error()
        assert max_error < 1e-10
        assert mae <= max_error
        rng = NumericalInversePolynomial(dist, u_resolution=1e-14)
        max_error, mae = rng.u_error()
        assert max_error < 1e-14
        assert mae <= max_error

    bad_orders = [1, 4.5, 20, np.inf, np.nan]
    bad_u_resolution = [1e-20, 1e-1, np.inf, np.nan]

    @pytest.mark.parametrize("order", bad_orders)
    def test_bad_orders(self, order):
        dist = StandardNormal()

        msg = r"`order` must be an integer in the range \[3, 17\]."
        with pytest.raises(ValueError, match=msg):
            NumericalInversePolynomial(dist, order=order)

    @pytest.mark.parametrize("u_resolution", bad_u_resolution)
    def test_bad_u_resolution(self, u_resolution):
        msg = r"`u_resolution` must be between 1e-15 and 1e-5."
        with pytest.raises(ValueError, match=msg):
            NumericalInversePolynomial(StandardNormal(),
                                       u_resolution=u_resolution)

    def test_bad_args(self):

        class BadDist:
            def cdf(self, x):
                return stats.norm._cdf(x)

        dist = BadDist()
        msg = r"Either of the methods `pdf` or `logpdf` must be specified"
        with pytest.raises(ValueError, match=msg):
            rng = NumericalInversePolynomial(dist)

        dist = StandardNormal()
        rng = NumericalInversePolynomial(dist)
        msg = r"`sample_size` must be greater than or equal to 1000."
        with pytest.raises(ValueError, match=msg):
            rng.u_error(10)

        class Distribution:
            def pdf(self, x):
                return np.exp(-0.5 * x*x)

        dist = Distribution()
        rng = NumericalInversePolynomial(dist)
        msg = r"Exact CDF required but not found."
        with pytest.raises(ValueError, match=msg):
            rng.u_error()

    def test_logpdf_pdf_consistency(self):
        # 1. check that PINV works with pdf and logpdf only
        # 2. check that generated ppf is the same (up to a small tolerance)

        class MyDist:
            pass

        # create generator from dist with only pdf
        dist_pdf = MyDist()
        dist_pdf.pdf = lambda x: math.exp(-x*x/2)
        rng1 = NumericalInversePolynomial(dist_pdf)

        # create dist with only logpdf
        dist_logpdf = MyDist()
        dist_logpdf.logpdf = lambda x: -x*x/2
        rng2 = NumericalInversePolynomial(dist_logpdf)

        q = np.linspace(1e-5, 1-1e-5, num=100)
        assert_allclose(rng1.ppf(q), rng2.ppf(q))


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])
    @pytest.mark.thread_unsafe
    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)

    def basic_test_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', 'invgauss', 'wald'}

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

        rng = np.random.default_rng(0)

        dist = getattr(stats, distname)(*shapes)
        fni = NumericalInverseHermite(dist)

        x = rng.random(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

    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
    @pytest.mark.xslow
    @pytest.mark.parametrize(("distname", "shapes"), distcont)
    def test_basic_all_scipy_dists(self, distname, shapes):
        # if distname == "truncnorm":
        #     pytest.skip("Tested separately")
        self.basic_test_all_scipy_dists(distname, shapes)

    @pytest.mark.fail_slow(5)
    @pytest.mark.filterwarnings('ignore::RuntimeWarning')
    @pytest.mark.parallel_threads_limit(4)  # Very slow
    def test_basic_truncnorm_gh17155(self):
        self.basic_test_all_scipy_dists("truncnorm", (0.1, 2))

    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):
            NumericalInverseHermite("norm")

        match = "could not convert string to float"
        with pytest.raises(ValueError, match=match):
            NumericalInverseHermite(StandardNormal(),
                                    u_resolution='ekki')

    rngs = [None, 0, np.random.RandomState(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)
    @pytest.mark.thread_unsafe
    def test_RVS(self, rng, size_in, size_out):
        dist = StandardNormal()
        fni = 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)

    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):
            NumericalInverseHermite(stats.beta(*shapes))

        # no error with coarser tol
        NumericalInverseHermite(stats.beta(*shapes), u_resolution=1e-8)

    def test_custom_distribution(self):
        dist1 = StandardNormal()
        fni1 = NumericalInverseHermite(dist1)

        dist2 = stats.norm()
        fni2 = 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 warnings.catch_warnings():
            msg = "invalid value encountered in "
            warnings.filterwarnings("ignore", msg + "greater", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "greater_equal", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less_equal", RuntimeWarning)
            res = rng.ppf(u)
            expected = stats.norm.ppf(u)
        assert_allclose(res, expected, rtol=1e-9, atol=3e-10)
        assert res.shape == expected.shape

    @pytest.mark.slow
    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 warnings.catch_warnings():
            # ignore warning about u-resolution being too small.
            warnings.simplefilter("ignore", RuntimeWarning)
            rng = NumericalInverseHermite(dist, u_resolution=1e-14)
        max_error, mae = rng.u_error()
        assert max_error < 1e-14
        assert mae <= max_error


class TestDiscreteGuideTable:
    basic_fail_dists = {
        'nchypergeom_fisher',  # numerical errors on tails
        'nchypergeom_wallenius',  # numerical errors on tails
        'randint'  # fails on 32-bit ubuntu
    }

    def test_guide_factor_gt3_raises_warning(self):
        pv = [0.1, 0.3, 0.6]
        urng = np.random.default_rng()
        with pytest.warns(RuntimeWarning):
            DiscreteGuideTable(pv, random_state=urng, guide_factor=7)

    def test_guide_factor_zero_raises_warning(self):
        pv = [0.1, 0.3, 0.6]
        urng = np.random.default_rng()
        with pytest.warns(RuntimeWarning):
            DiscreteGuideTable(pv, random_state=urng, guide_factor=0)

    def test_negative_guide_factor_raises_warning(self):
        # This occurs from the UNU.RAN wrapper automatically.
        # however it already gives a useful warning
        # Here we just test that a warning is raised.
        pv = [0.1, 0.3, 0.6]
        urng = np.random.default_rng()
        with pytest.warns(RuntimeWarning):
            DiscreteGuideTable(pv, random_state=urng, guide_factor=-1)

    @pytest.mark.parametrize("distname, params", distdiscrete)
    def test_basic(self, distname, params):
        if distname in self.basic_fail_dists:
            msg = ("DGT fails on these probably because of large domains "
                   "and small computation errors in PMF.")
            pytest.skip(msg)

        if not isinstance(distname, str):
            dist = distname
        else:
            dist = getattr(stats, distname)

        dist = dist(*params)
        domain = dist.support()

        if not np.isfinite(domain[1] - domain[0]):
            # DGT only works with finite domain. So, skip the distributions
            # with infinite tails.
            pytest.skip("DGT only works with a finite domain.")

        k = np.arange(domain[0], domain[1]+1)
        pv = dist.pmf(k)
        mv_ex = dist.stats('mv')
        rng = DiscreteGuideTable(dist, random_state=42)
        check_discr_samples(rng, pv, mv_ex)

    u = [
        # the correctness of the PPF for equidistant points between 0 and 1.
        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):
        n, p = 4, 0.1
        dist = stats.binom(n, p)
        rng = DiscreteGuideTable(dist, random_state=42)

        # Older versions of NumPy throw RuntimeWarnings for comparisons
        # with nan.
        with warnings.catch_warnings():
            msg = "invalid value encountered in "
            warnings.filterwarnings("ignore", msg + "greater", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "greater_equal", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less", RuntimeWarning)
            warnings.filterwarnings("ignore", msg + "less_equal", RuntimeWarning)
            res = rng.ppf(u)
            expected = stats.binom.ppf(u, n, p)
        assert_equal(res.shape, expected.shape)
        assert_equal(res, expected)

    @pytest.mark.parametrize("pv, msg", bad_pv_common)
    def test_bad_pv(self, pv, msg):
        with pytest.raises(ValueError, match=msg):
            DiscreteGuideTable(pv)

    # DGT doesn't support infinite tails. So, it should throw an error when
    # inf is present in the domain.
    inf_domain = [(-np.inf, np.inf), (np.inf, np.inf), (-np.inf, -np.inf),
                  (0, np.inf), (-np.inf, 0)]

    @pytest.mark.parametrize("domain", inf_domain)
    def test_inf_domain(self, domain):
        with pytest.raises(ValueError, match=r"must be finite"):
            DiscreteGuideTable(stats.binom(10, 0.2), domain=domain)


class TestSimpleRatioUniforms:
    # pdf with piecewise linear function as transformed density
    # with T = -1/sqrt with shift. Taken from UNU.RAN test suite
    # (from file t_srou.c)
    class dist:
        def __init__(self, shift):
            self.shift = shift
            self.mode = shift

        def pdf(self, x):
            x -= self.shift
            y = 1. / (abs(x) + 1.)
            return 0.5 * y * y

        def cdf(self, x):
            x -= self.shift
            if x <= 0.:
                return 0.5 / (1. - x)
            else:
                return 1. - 0.5 / (1. + x)

    dists = [dist(0.), dist(10000.)]

    # exact mean and variance of the distributions in the list dists
    mv1 = [0., np.inf]
    mv2 = [10000., np.inf]
    mvs = [mv1, mv2]

    @pytest.mark.parametrize("dist, mv_ex",
                             zip(dists, mvs))
    @pytest.mark.thread_unsafe
    def test_basic(self, dist, mv_ex):
        rng = SimpleRatioUniforms(dist, mode=dist.mode, random_state=42)
        check_cont_samples(rng, dist, mv_ex)
        rng = SimpleRatioUniforms(dist, mode=dist.mode,
                                  cdf_at_mode=dist.cdf(dist.mode),
                                  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):
            SimpleRatioUniforms(StandardNormal(), domain=domain)

    def test_bad_args(self):
        # pdf_area < 0
        with pytest.raises(ValueError, match=r"`pdf_area` must be > 0"):
            SimpleRatioUniforms(StandardNormal(), mode=0, pdf_area=-1)


class TestRatioUniforms:
    def test_rv_generation(self):
        # use KS test to check distribution of rvs
        # normal distribution
        f = stats.norm.pdf
        v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
        u = np.sqrt(f(0))
        gen = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=12345)
        assert_equal(stats.kstest(gen.rvs(2500), 'norm')[1] > 0.25, True)

        # exponential distribution
        gen = RatioUniforms(lambda x: np.exp(-x), umax=1,
                            vmin=0, vmax=2*np.exp(-1), random_state=12345)
        assert_equal(stats.kstest(gen.rvs(1000), 'expon')[1] > 0.25, True)

    def test_shape(self):
        # test shape of return value depending on size parameter
        f = stats.norm.pdf
        v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
        u = np.sqrt(f(0))

        gen1 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        gen2 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        gen3 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        r1, r2, r3 = gen1.rvs(3), gen2.rvs((3,)), gen3.rvs((3, 1))
        assert_equal(r1, r2)
        assert_equal(r2, r3.flatten())
        assert_equal(r1.shape, (3,))
        assert_equal(r3.shape, (3, 1))

        gen4 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=12)
        gen5 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=12)
        r4, r5 = gen4.rvs(size=(3, 3, 3)), gen5.rvs(size=27)
        assert_equal(r4.flatten(), r5)
        assert_equal(r4.shape, (3, 3, 3))

        gen6 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        gen7 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        gen8 = RatioUniforms(f, umax=u, vmin=-v, vmax=v, random_state=1234)
        r6, r7, r8 = gen6.rvs(), gen7.rvs(1), gen8.rvs((1,))
        assert_equal(r6, r7)
        assert_equal(r7, r8)

    def test_random_state(self):
        f = stats.norm.pdf
        v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
        umax = np.sqrt(f(0))
        gen1 = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=1234)
        r1 = gen1.rvs(10)
        rng = np.random.RandomState(1234)
        gen2 = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng)
        r2 = gen2.rvs(10)
        assert_equal(r1, r2)

    def test_exceptions(self):
        f = stats.norm.pdf
        # need vmin < vmax
        with assert_raises(ValueError, match="vmin must be smaller than vmax"):
            RatioUniforms(pdf=f, umax=1, vmin=3, vmax=1)
        with assert_raises(ValueError, match="vmin must be smaller than vmax"):
            RatioUniforms(pdf=f, umax=1, vmin=1, vmax=1)
        # need umax > 0
        with assert_raises(ValueError, match="umax must be positive"):
            RatioUniforms(pdf=f, umax=-1, vmin=1, vmax=3)
        with assert_raises(ValueError, match="umax must be positive"):
            RatioUniforms(pdf=f, umax=0, vmin=1, vmax=3)