diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 6111ae464a..398ee38b84 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -9,6 +9,7 @@ - The `Distribution` keyword argument `testval` has been deprecated in favor of `initval`. Furthermore `initval` no longer assigns a `tag.test_value` on tensors since the initial values are now kept track of by the model object ([see #4913](https://github.com/pymc-devs/pymc3/pull/4913)). - `pm.sample` now returns results as `InferenceData` instead of `MultiTrace` by default (see [#4744](https://github.com/pymc-devs/pymc3/pull/4744)). - `pm.sample_prior_predictive` no longer returns transformed variable values by default. Pass them by name in `var_names` if you want to obtain these draws (see [4769](https://github.com/pymc-devs/pymc3/pull/4769)). +- ⚠ `pm.Bound` interface no longer accepts a callable class as argument, instead it requires an instantiated distribution (created via the `.dist()` API) to be passed as an argument. In addition, Bound no longer returns a class instance but works as a normal PyMC3 distribution. Finally, it is no longer possible to do predictive random sampling from Bounded variables. Please, consult the new documentation for details on how to use Bounded variables (see [4815](https://github.com/pymc-devs/pymc3/pull/4815)). - ... ### New Features diff --git a/docs/source/api/bounds.rst b/docs/source/api/bounds.rst index 250044cfe9..d4a0b003a6 100644 --- a/docs/source/api/bounds.rst +++ b/docs/source/api/bounds.rst @@ -6,10 +6,6 @@ PyMC3 includes the construct :class:`~pymc3.distributions.bound.Bound` for placing constraints on existing probability distributions. It modifies a given distribution to take values only within a specified interval. -Note that the `Bound` class does *not* directly create a bounded -distribution: instead it creates a Callable class that can be -*invoked* to create a bounded distribution, as the example below illustrates. - Some types of variables require constraints. For instance, it doesn't make sense for a standard deviation to have a negative value, so something like a Normal prior on a parameter that represents a standard deviation would be @@ -38,28 +34,8 @@ specification of a bounded distribution should go within the model block:: import pymc3 as pm with pm.Model() as model: - BoundedNormal = pm.Bound(pm.Normal, lower=0.0) - x = BoundedNormal('x', mu=1.0, sigma=3.0) - -If the bound will be applied to a single variable in the model, it may be -cleaner notationally to define both the bound and variable together. :: - - with model: - x = pm.Bound(pm.Normal, lower=0.0)('x', mu=1.0, sigma=3.0) - -However, it is possible to create multiple different random variables -that have the same bound applied to them:: - - with model: - BoundNormal = pm.Bound(pm.Normal, lower=0.0) - hyper_mu = BoundNormal("hyper_mu", mu=1, sigma=0.5) - mu = BoundNormal("mu", mu=hyper_mu, sigma=1) - -Bounds can also be applied to a vector of random variables. With the same -``BoundedNormal`` object we created previously we can write:: - - with model: - x_vector = BoundedNormal('x_vector', mu=1.0, sigma=3.0, shape=3) + norm = Normal.dist(mu=1.0, sigma=3.0) + x = Bound('x', norm, lower=0.0) Caveats ####### diff --git a/pymc3/distributions/bound.py b/pymc3/distributions/bound.py index 89f02bdcfb..19169d75b5 100644 --- a/pymc3/distributions/bound.py +++ b/pymc3/distributions/bound.py @@ -12,42 +12,43 @@ # See the License for the specific language governing permissions and # limitations under the License. -from numbers import Real - -import aesara.tensor as at import numpy as np -from pymc3.aesaraf import floatX -from pymc3.distributions import transforms +from aesara.tensor import as_tensor_variable +from aesara.tensor.random.op import RandomVariable +from aesara.tensor.var import TensorVariable + +from pymc3.aesaraf import floatX, intX +from pymc3.distributions import _logp +from pymc3.distributions.continuous import BoundedContinuous from pymc3.distributions.dist_math import bound -from pymc3.distributions.distribution import Continuous, Discrete, Distribution +from pymc3.distributions.distribution import Continuous, Discrete +from pymc3.distributions.shape_utils import to_tuple +from pymc3.model import modelcontext __all__ = ["Bound"] -class _Bounded(Distribution): - def __init__(self, distribution, lower, upper, default, *args, **kwargs): - self.lower = lower - self.upper = upper - self._wrapped = distribution.dist(*args, **kwargs) +class BoundRV(RandomVariable): + name = "bound" + ndim_supp = 0 + ndims_params = [0, 0, 0] + dtype = "floatX" + _print_name = ("Bound", "\\operatorname{Bound}") - if default is None: - defaults = self._wrapped.defaults - for name in defaults: - setattr(self, name, getattr(self._wrapped, name)) - else: - defaults = ("_default",) - self._default = default - - super().__init__( - shape=self._wrapped.shape, - dtype=self._wrapped.dtype, - initval=self._wrapped.initval, - defaults=defaults, - transform=self._wrapped.transform, - ) - - def logp(self, value): + @classmethod + def rng_fn(cls, rng, distribution, lower, upper, size): + raise NotImplementedError("Cannot sample from a bounded variable") + + +boundrv = BoundRV() + + +class _ContinuousBounded(BoundedContinuous): + rv_op = boundrv + bound_args_indices = [1, 2] + + def logp(value, distribution, lower, upper): """ Calculate log-probability of Bounded distribution at specified value. @@ -55,155 +56,59 @@ def logp(self, value): ---------- value: numeric Value for which log-probability is calculated. + distribution: TensorVariable + Distribution which is being bounded + lower: numeric + Lower bound for the distribution being bounded. + upper: numeric + Upper bound for the distribution being bounded. Returns ------- TensorVariable """ - logp = self._wrapped.logp(value) - bounds = [] - if self.lower is not None: - bounds.append(value >= self.lower) - if self.upper is not None: - bounds.append(value <= self.upper) - if len(bounds) > 0: - return bound(logp, *bounds) - else: - return logp - - def _random(self, lower, upper, point=None, size=None): - lower = np.asarray(lower) - upper = np.asarray(upper) - if lower.size > 1 or upper.size > 1: - raise ValueError( - "Drawing samples from distributions with array-valued bounds is not supported." - ) - total_size = np.prod(size).astype(np.int) - samples = [] - s = 0 - while s < total_size: - sample = np.atleast_1d(self._wrapped.random(point=point, size=total_size)).flatten() - - select = sample[np.logical_and(sample >= lower, sample <= upper)] - samples.append(select) - s += len(select) - if size is not None: - return np.reshape(np.concatenate(samples)[:total_size], size) - else: - return samples[0] + logp = _logp(distribution.owner.op, value, {}, *distribution.owner.inputs[3:]) + return bound(logp, (value >= lower), (value <= upper)) - def random(self, point=None, size=None): - """ - Draw random values from Bounded distribution. - Parameters - ---------- - point: dict, optional - Dict of variable values on which random values are to be - conditioned (uses default point if not specified). - size: int, optional - Desired size of random sample (returns one sample if not - specified). +class DiscreteBoundRV(BoundRV): + name = "discrete_bound" + dtype = "int64" - Returns - ------- - array - """ - # if self.lower is None and self.upper is None: - # return self._wrapped.random(point=point, size=size) - # elif self.lower is not None and self.upper is not None: - # lower, upper = draw_values([self.lower, self.upper], point=point, size=size) - # return generate_samples( - # self._random, - # lower, - # upper, - # dist_shape=self.shape, - # size=size, - # not_broadcast_kwargs={"point": point}, - # ) - # elif self.lower is not None: - # lower = draw_values([self.lower], point=point, size=size) - # return generate_samples( - # self._random, - # lower, - # np.inf, - # dist_shape=self.shape, - # size=size, - # not_broadcast_kwargs={"point": point}, - # ) - # else: - # upper = draw_values([self.upper], point=point, size=size) - # return generate_samples( - # self._random, - # -np.inf, - # upper, - # dist_shape=self.shape, - # size=size, - # not_broadcast_kwargs={"point": point}, - # ) - pass - - -class _DiscreteBounded(_Bounded, Discrete): - def __init__(self, distribution, lower, upper, transform="infer", *args, **kwargs): - if transform == "infer": - transform = None - if transform is not None: - raise ValueError("Can not transform discrete variable.") - if lower is None and upper is None: - default = None - elif lower is not None and upper is not None: - default = (lower + upper) // 2 - if upper is not None: - default = upper - 1 - if lower is not None: - default = lower + 1 +discrete_boundrv = DiscreteBoundRV() - super().__init__(distribution, lower, upper, default, *args, transform=transform, **kwargs) +class _DiscreteBounded(Discrete): + rv_op = discrete_boundrv -class _ContinuousBounded(_Bounded, Continuous): - r""" - An upper, lower or upper+lower bounded distribution + def __new__(cls, *args, **kwargs): + transform = kwargs.get("transform", None) + if transform is not None: + raise ValueError("Cannot transform discrete variable.") + return super().__new__(cls, *args, **kwargs) - Parameters - ---------- - distribution: pymc3 distribution - Distribution to be transformed into a bounded distribution - lower: float (optional) - Lower bound of the distribution, set to -inf to disable. - upper: float (optional) - Upper bound of the distribibution, set to inf to disable. - tranform: 'infer' or object - If 'infer', infers the right transform to apply from the supplied bounds. - If transform object, has to supply .forward() and .backward() methods. - See pymc3.distributions.transforms for more information. - """ + def logp(value, distribution, lower, upper): + """ + Calculate log-probability of Bounded distribution at specified value. - def __init__(self, distribution, lower, upper, transform="infer", *args, **kwargs): - if lower is not None: - lower = at.as_tensor_variable(floatX(lower)) - if upper is not None: - upper = at.as_tensor_variable(floatX(upper)) - - if transform == "infer": - if lower is None and upper is None: - transform = None - default = None - elif lower is not None and upper is not None: - transform = transforms.interval(lower, upper) - default = 0.5 * (lower + upper) - elif upper is not None: - transform = transforms.upperbound(upper) - default = upper - 1 - else: - transform = transforms.lowerbound(lower) - default = lower + 1 - else: - default = None + Parameters + ---------- + value: numeric + Value for which log-probability is calculated. + distribution: TensorVariable + Distribution which is being bounded + lower: numeric + Lower bound for the distribution being bounded. + upper: numeric + Upper bound for the distribution being bounded. - super().__init__(distribution, lower, upper, default, *args, transform=transform, **kwargs) + Returns + ------- + TensorVariable + """ + logp = _logp(distribution.owner.op, value, {}, *distribution.owner.inputs[3:]) + return bound(logp, (value >= lower), (value <= upper)) class Bound: @@ -232,32 +137,90 @@ class Bound: .. code-block:: python with pm.Model(): - NegativeNormal = pm.Bound(pm.Normal, upper=0.0) - par1 = NegativeNormal('par`', mu=0.0, sigma=1.0, initval=-0.5) - # you can use the Bound object multiple times to - # create multiple bounded random variables - par1_1 = NegativeNormal('par1_1', mu=-1.0, sigma=1.0, initval=-1.5) - - # you can also define a Bound implicitly, while applying - # it to a random variable - par2 = pm.Bound(pm.Normal, lower=-1.0, upper=1.0)( - 'par2', mu=0.0, sigma=1.0, initval=1.0) - """ - - def __init__(self, distribution, lower=None, upper=None): - if isinstance(lower, Real) and lower == -np.inf: - lower = None - if isinstance(upper, Real) and upper == np.inf: - upper = None + normal_dist = Normal.dist(mu=0.0, sigma=1.0, initval=-0.5) + negative_normal = pm.Bound(normal_dist, upper=0.0) - if not issubclass(distribution, Distribution): - raise ValueError('"distribution" must be a Distribution subclass.') + """ - self.distribution = distribution - self.lower = lower - self.upper = upper + def __new__( + cls, + name, + distribution, + lower=None, + upper=None, + size=None, + shape=None, + initval=None, + dims=None, + **kwargs, + ): + + cls._argument_checks(distribution, **kwargs) + + if dims is not None: + model = modelcontext(None) + if dims in model.coords: + dim_obj = np.asarray(model.coords[dims]) + size = dim_obj.shape + else: + raise ValueError("Given dims do not exist in model coordinates.") + + lower, upper, initval = cls._set_values(lower, upper, size, shape, initval) + + if isinstance(distribution.owner.op, Continuous): + res = _ContinuousBounded( + name, + [distribution, lower, upper], + initval=floatX(initval), + size=size, + shape=shape, + **kwargs, + ) + else: + res = _DiscreteBounded( + name, + [distribution, lower, upper], + initval=intX(initval), + size=size, + shape=shape, + **kwargs, + ) + return res + + @classmethod + def dist( + cls, + distribution, + lower=None, + upper=None, + size=None, + shape=None, + **kwargs, + ): + + cls._argument_checks(distribution, **kwargs) + lower, upper, initval = cls._set_values(lower, upper, size, shape, initval=None) + + if isinstance(distribution.owner.op, Continuous): + res = _ContinuousBounded.dist( + [distribution, lower, upper], + size=size, + shape=shape, + **kwargs, + ) + res.tag.test_value = floatX(initval) + else: + res = _DiscreteBounded.dist( + [distribution, lower, upper], + size=size, + shape=shape, + **kwargs, + ) + res.tag.test_value = intX(initval) + return res - def __call__(self, name, *args, **kwargs): + @classmethod + def _argument_checks(cls, distribution, **kwargs): if "observed" in kwargs: raise ValueError( "Observed Bound distributions are not supported. " @@ -266,25 +229,60 @@ def __call__(self, name, *args, **kwargs): "with the cumulative probability function." ) - transform = kwargs.pop("transform", "infer") - if issubclass(self.distribution, Continuous): - return _ContinuousBounded( - name, self.distribution, self.lower, self.upper, transform, *args, **kwargs - ) - elif issubclass(self.distribution, Discrete): - return _DiscreteBounded( - name, self.distribution, self.lower, self.upper, transform, *args, **kwargs + if not isinstance(distribution, TensorVariable): + raise ValueError( + "Passing a distribution class to `Bound` is no longer supported.\n" + "Please pass the output of a distribution instantiated via the " + "`.dist()` API such as:\n" + '`pm.Bound("bound", pm.Normal.dist(0, 1), lower=0)`' ) + + try: + model = modelcontext(None) + except TypeError: + pass else: - raise ValueError("Distribution is neither continuous nor discrete.") + if distribution in model.basic_RVs: + raise ValueError( + f"The distribution passed into `Bound` was already registered " + f"in the current model.\nYou should pass an unregistered " + f"(unnamed) distribution created via the `.dist()` API, such as:\n" + f'`pm.Bound("bound", pm.Normal.dist(0, 1), lower=0)`' + ) + + if distribution.owner.op.ndim_supp != 0: + raise NotImplementedError("Bounding of MultiVariate RVs is not yet supported.") + + if not isinstance(distribution.owner.op, (Discrete, Continuous)): + raise ValueError( + f"`distribution` {distribution} must be a Discrete or Continuous" + " distribution subclass" + ) + + @classmethod + def _set_values(cls, lower, upper, size, shape, initval): + if size is None: + size = shape - def dist(self, *args, **kwargs): - if issubclass(self.distribution, Continuous): - return _ContinuousBounded.dist( - self.distribution, self.lower, self.upper, *args, **kwargs + lower = np.asarray(lower) + lower = floatX(np.where(lower == None, -np.inf, lower)) + upper = np.asarray(upper) + upper = floatX(np.where(upper == None, np.inf, upper)) + + if initval is None: + _size = np.broadcast_shapes(to_tuple(size), np.shape(lower), np.shape(upper)) + _lower = np.broadcast_to(lower, _size) + _upper = np.broadcast_to(upper, _size) + initval = np.where( + (_lower == -np.inf) & (_upper == np.inf), + 0, + np.where( + _lower == -np.inf, + _upper - 1, + np.where(_upper == np.inf, _lower + 1, (_lower + _upper) / 2), + ), ) - elif issubclass(self.distribution, Discrete): - return _DiscreteBounded.dist(self.distribution, self.lower, self.upper, *args, **kwargs) - else: - raise ValueError("Distribution is neither continuous nor discrete.") + lower = as_tensor_variable(floatX(lower)) + upper = as_tensor_variable(floatX(upper)) + return lower, upper, initval diff --git a/pymc3/model.py b/pymc3/model.py index df3c486338..4b7df58c09 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -963,7 +963,7 @@ def set_initval(self, rv_var, initval): if transform: value = initval if initval is not None else rv_var - rv_var = transform.forward(rv_var, value) + rv_var = at.as_tensor_variable(transform.forward(rv_var, value)) def initval_to_rvval(value_var, value): rv_var = self.values_to_rvs[value_var] diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 261aeae2bd..5c46e14105 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -69,6 +69,7 @@ def polyagamma_cdf(*args, **kwargs): DirichletMultinomial, DiscreteUniform, DiscreteWeibull, + Distribution, ExGaussian, Exponential, Flat, @@ -956,20 +957,6 @@ def test_triangular(self): assert logp(invalid_dist, 0.5).eval() == -np.inf assert logcdf(invalid_dist, 2).eval() == -np.inf - @pytest.mark.xfail(reason="Bound not refactored yet") - def test_bound_normal(self): - PositiveNormal = Bound(Normal, lower=0.0) - self.check_logp( - PositiveNormal, - Rplus, - {"mu": Rplus, "sigma": Rplus}, - lambda value, mu, sigma: sp.norm.logpdf(value, mu, sigma), - decimal=select_by_precision(float64=6, float32=-1), - ) - with Model(): - x = PositiveNormal("x", mu=0, sigma=1, transform=None) - assert np.isinf(logp(x, -1).eval()) - @pytest.mark.skipif( condition=_polyagamma_not_installed, reason="`polyagamma package is not available/installed.", @@ -1687,20 +1674,6 @@ def test_poisson(self): {"mu": Rplus}, ) - @pytest.mark.xfail(reason="Distribution not refactored yet") - def test_bound_poisson(self): - NonZeroPoisson = Bound(Poisson, lower=1.0) - self.check_logp( - NonZeroPoisson, - PosNat, - {"mu": Rplus}, - lambda value, mu: sp.poisson.logpmf(value, mu), - ) - - with Model(): - x = NonZeroPoisson("x", mu=4) - assert np.isinf(logp(x, 0).eval()) - def test_constantdist(self): self.check_logp(Constant, I, {"c": I}, lambda value, c: np.log(c == value)) @@ -2669,74 +2642,184 @@ def ref_pdf(value): self.check_logp(TestedInterpolated, R, {}, ref_pdf) -@pytest.mark.xfail(reason="Bound not refactored yet") -def test_bound(): - np.random.seed(42) - UnboundNormal = Bound(Normal) - dist = UnboundNormal.dist(mu=0, sigma=1) - # assert dist.transform is None - assert isinstance(dist.random(), np.ndarray) - - LowerNormal = Bound(Normal, lower=1) - dist = LowerNormal.dist(mu=0, sigma=1) - assert logp(dist, 0).eval() == -np.inf - # assert dist.transform is not None - assert np.all(dist.random() > 1) - - UpperNormal = Bound(Normal, upper=-1) - dist = UpperNormal.dist(mu=0, sigma=1) - assert logp(dist, -0.5).eval() == -np.inf - # assert dist.transform is not None - assert np.all(dist.random() < -1) - - ArrayNormal = Bound(Normal, lower=[1, 2], upper=[2, 3]) - dist = ArrayNormal.dist(mu=0, sigma=1, size=2) - assert_equal(logp(dist, [0.5, 3.5]).eval(), -np.array([np.inf, np.inf])) - # assert dist.transform is not None - with pytest.raises(ValueError) as err: - dist.random() - err.match("Drawing samples from distributions with array-valued") - - with Model(): - a = ArrayNormal("c", size=2) - assert_equal(a.tag.test_value, np.array([1.5, 2.5])) - - lower = at.vector("lower") - lower.tag.test_value = np.array([1, 2]).astype(aesara.config.floatX) - upper = 3 - ArrayNormal = Bound(Normal, lower=lower, upper=upper) - dist = ArrayNormal.dist(mu=0, sigma=1, size=2) - logp_dist = logp(dist, [0.5, 3.5]).eval({lower: lower.tag.test_value}) - assert_equal(logp_dist, -np.array([np.inf, np.inf])) - assert dist.transform is not None - - with Model(): - a = ArrayNormal("c", size=2) - assert_equal(a.tag.test_value, np.array([2, 2.5])) - - rand = Bound(Binomial, lower=10).dist(n=20, p=0.3).random() - assert rand.dtype in [np.int16, np.int32, np.int64] - assert rand >= 10 - - rand = Bound(Binomial, upper=10).dist(n=20, p=0.8).random() - assert rand.dtype in [np.int16, np.int32, np.int64] - assert rand <= 10 - - rand = Bound(Binomial, lower=5, upper=8).dist(n=10, p=0.6).random() - assert rand.dtype in [np.int16, np.int32, np.int64] - assert rand >= 5 and rand <= 8 +class TestBound: + """Tests for pm.Bound distribution""" - with Model(): - BoundPoisson = Bound(Poisson, upper=6) - BoundPoisson(name="y", mu=1) - - with Model(): - BoundNormalNamedArgs = Bound(Normal, upper=6)("y", mu=2.0, sd=1.0) - BoundNormalPositionalArgs = Bound(Normal, upper=6)("x", 2.0, 1.0) - - with Model(): - BoundPoissonNamedArgs = Bound(Poisson, upper=6)("y", mu=2.0) - BoundPoissonPositionalArgs = Bound(Poisson, upper=6)("x", 2.0) + def test_continuous(self): + with Model() as model: + dist = Normal.dist(mu=0, sigma=1) + UnboundedNormal = Bound("unbound", dist, transform=None) + InfBoundedNormal = Bound("infbound", dist, lower=-np.inf, upper=np.inf, transform=None) + LowerNormal = Bound("lower", dist, lower=0, transform=None) + UpperNormal = Bound("upper", dist, upper=0, transform=None) + BoundedNormal = Bound("bounded", dist, lower=1, upper=10, transform=None) + LowerNormalTransform = Bound("lowertrans", dist, lower=1) + UpperNormalTransform = Bound("uppertrans", dist, upper=10) + BoundedNormalTransform = Bound("boundedtrans", dist, lower=1, upper=10) + + assert logpt(LowerNormal, -1).eval() == -np.inf + assert logpt(UpperNormal, 1).eval() == -np.inf + assert logpt(BoundedNormal, 0).eval() == -np.inf + assert logpt(BoundedNormal, 11).eval() == -np.inf + + assert logpt(UnboundedNormal, 0).eval() != -np.inf + assert logpt(UnboundedNormal, 11).eval() != -np.inf + assert logpt(InfBoundedNormal, 0).eval() != -np.inf + assert logpt(InfBoundedNormal, 11).eval() != -np.inf + + assert logpt(LowerNormalTransform, -1).eval() != -np.inf + assert logpt(UpperNormalTransform, 1).eval() != -np.inf + assert logpt(BoundedNormalTransform, 0).eval() != -np.inf + assert logpt(BoundedNormalTransform, 11).eval() != -np.inf + + assert np.allclose( + logpt(UnboundedNormal, 5).eval(), Normal.logp(value=5, mu=0, sigma=1).eval() + ) + assert np.allclose(logpt(LowerNormal, 5).eval(), Normal.logp(value=5, mu=0, sigma=1).eval()) + assert np.allclose( + logpt(UpperNormal, -5).eval(), Normal.logp(value=-5, mu=0, sigma=1).eval() + ) + assert np.allclose( + logpt(BoundedNormal, 5).eval(), Normal.logp(value=5, mu=0, sigma=1).eval() + ) + + def test_discrete(self): + with Model() as model: + dist = Poisson.dist(mu=4) + UnboundedPoisson = Bound("unbound", dist) + LowerPoisson = Bound("lower", dist, lower=1) + UpperPoisson = Bound("upper", dist, upper=10) + BoundedPoisson = Bound("bounded", dist, lower=1, upper=10) + + assert logpt(LowerPoisson, 0).eval() == -np.inf + assert logpt(UpperPoisson, 11).eval() == -np.inf + assert logpt(BoundedPoisson, 0).eval() == -np.inf + assert logpt(BoundedPoisson, 11).eval() == -np.inf + + assert logpt(UnboundedPoisson, 0).eval() != -np.inf + assert logpt(UnboundedPoisson, 11).eval() != -np.inf + + assert np.allclose(logpt(UnboundedPoisson, 5).eval(), Poisson.logp(value=5, mu=4).eval()) + assert np.allclose(logpt(UnboundedPoisson, 5).eval(), Poisson.logp(value=5, mu=4).eval()) + assert np.allclose(logpt(UnboundedPoisson, 5).eval(), Poisson.logp(value=5, mu=4).eval()) + assert np.allclose(logpt(UnboundedPoisson, 5).eval(), Poisson.logp(value=5, mu=4).eval()) + + def create_invalid_distribution(self): + class MyNormal(RandomVariable): + name = "my_normal" + ndim_supp = 0 + ndims_params = [0, 0] + dtype = "floatX" + + my_normal = MyNormal() + + class InvalidDistribution(Distribution): + rv_op = my_normal + + @classmethod + def dist(cls, mu=0, sigma=1, **kwargs): + return super().dist([mu, sigma], **kwargs) + + return InvalidDistribution + + def test_arguments_checks(self): + msg = "Observed Bound distributions are not supported" + with pm.Model() as m: + x = pm.Normal("x", 0, 1) + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", x, observed=5) + + msg = "Cannot transform discrete variable." + with pm.Model() as m: + x = pm.Poisson.dist(0.5) + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", x, transform=pm.transforms.interval) + + msg = "Given dims do not exist in model coordinates." + with pm.Model() as m: + x = pm.Poisson.dist(0.5) + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", x, dims="random_dims") + + msg = "The distribution passed into `Bound` was already registered" + with pm.Model() as m: + x = pm.Normal("x", 0, 1) + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", x) + + msg = "Passing a distribution class to `Bound` is no longer supported" + with pm.Model() as m: + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", pm.Normal) + + msg = "Bounding of MultiVariate RVs is not yet supported" + with pm.Model() as m: + x = pm.MvNormal.dist(np.zeros(3), np.eye(3)) + with pytest.raises(NotImplementedError, match=msg): + pm.Bound("bound", x) + + msg = "must be a Discrete or Continuous distribution subclass" + with pm.Model() as m: + x = self.create_invalid_distribution().dist() + with pytest.raises(ValueError, match=msg): + pm.Bound("bound", x) + + def test_invalid_sampling(self): + msg = "Cannot sample from a bounded variable" + with pm.Model() as m: + dist = Normal.dist(mu=0, sigma=1) + BoundedNormal = Bound("bounded", dist, lower=1, upper=10) + with pytest.raises(NotImplementedError, match=msg): + pm.sample_prior_predictive() + + def test_bound_shapes(self): + with pm.Model(coords={"sample": np.ones((2, 5))}) as m: + dist = Normal.dist(mu=0, sigma=1) + bound_sized = Bound("boundedsized", dist, lower=1, upper=10, size=(4, 5)) + bound_shaped = Bound("boundedshaped", dist, lower=1, upper=10, shape=(3, 5)) + bound_dims = Bound("boundeddims", dist, lower=1, upper=10, dims="sample") + + dist_size = m.initial_point["boundedsized_interval__"].shape + dist_shape = m.initial_point["boundedshaped_interval__"].shape + dist_dims = m.initial_point["boundeddims_interval__"].shape + + assert dist_size == (4, 5) + assert dist_shape == (3, 5) + assert dist_dims == (2, 5) + + def test_bound_dist(self): + # Continuous + bound = pm.Bound.dist(pm.Normal.dist(0, 1), lower=0) + assert pm.logp(bound, -1).eval() == -np.inf + assert np.isclose(pm.logp(bound, 1).eval(), scipy.stats.norm(0, 1).logpdf(1)) + + # Discrete + bound = pm.Bound.dist(pm.Poisson.dist(1), lower=2) + assert pm.logp(bound, 1).eval() == -np.inf + assert np.isclose(pm.logp(bound, 2).eval(), scipy.stats.poisson(1).logpmf(2)) + + def test_array_bound(self): + with Model() as model: + dist = Normal.dist() + LowerPoisson = Bound("lower", dist, lower=[1, None], transform=None) + UpperPoisson = Bound("upper", dist, upper=[np.inf, 10], transform=None) + BoundedPoisson = Bound("bounded", dist, lower=[1, 2], upper=[9, 10], transform=None) + + first, second = logpt(LowerPoisson, [0, 0]).eval() + assert first == -np.inf + assert second != -np.inf + + first, second = logpt(UpperPoisson, [11, 11]).eval() + assert first != -np.inf + assert second == -np.inf + + first, second = logpt(BoundedPoisson, [1, 1]).eval() + assert first != -np.inf + assert second == -np.inf + + first, second = logpt(BoundedPoisson, [10, 10]).eval() + assert first == -np.inf + assert second != -np.inf class TestBoundedContinuous: diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 9e6af20996..3c559b775c 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -1677,16 +1677,6 @@ def kronecker_rng_fn(self, size, mu, covs=None, sigma=None, rng=None): class TestScalarParameterSamples(SeededTest): - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") - def test_bounded(self): - # A bit crude... - BoundedNormal = pm.Bound(pm.Normal, upper=0) - - def ref_rand(size, tau): - return -st.halfnorm.rvs(size=size, loc=0, scale=tau ** -0.5) - - pymc3_random(BoundedNormal, {"tau": Rplus}, ref_rand=ref_rand) - @pytest.mark.xfail(reason="This distribution has not been refactored for v4") def test_lkj(self): for n in [2, 10, 50]: diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py index 80335bf7a4..01e28c9155 100644 --- a/pymc3/tests/test_sampling.py +++ b/pymc3/tests/test_sampling.py @@ -1058,16 +1058,6 @@ def test_zeroinflatedpoisson(self): assert gen_data["psi"].shape == (5000,) assert gen_data["suppliers"].shape == (5000, 20) - @pytest.mark.xfail(reason="Bound not refactored for v4") - def test_bounded_dist(self): - with pm.Model() as model: - BoundedNormal = pm.Bound(pm.Normal, lower=0.0) - x = BoundedNormal("x", mu=at.zeros((3, 1)), sd=1 * at.ones((3, 1)), size=(3, 1)) - - with model: - prior_trace = pm.sample_prior_predictive(5) - assert prior_trace["x"].shape == (5, 3, 1) - def test_potentials_warning(self): warning_msg = "The effect of Potentials on other parameters is ignored during" with pm.Model() as m: