diff --git a/pymc3/distributions/__init__.py b/pymc3/distributions/__init__.py index f8390db75f..65f0664f6d 100644 --- a/pymc3/distributions/__init__.py +++ b/pymc3/distributions/__init__.py @@ -42,7 +42,6 @@ from .discrete import Geometric from .discrete import Categorical -from .distribution import Bound from .distribution import DensityDist from .distribution import Distribution from .distribution import Continuous @@ -76,6 +75,8 @@ from .transforms import log from .transforms import sum_to_1 +from .bound import Bound + __all__ = ['Uniform', 'Flat', 'Normal', @@ -136,5 +137,6 @@ 'Triangular', 'DiscreteWeibull', 'Gumbel', - 'Interpolated' + 'Interpolated', + 'Bound', ] diff --git a/pymc3/distributions/bound.py b/pymc3/distributions/bound.py new file mode 100644 index 0000000000..8dc4e7966f --- /dev/null +++ b/pymc3/distributions/bound.py @@ -0,0 +1,232 @@ +from numbers import Real + +import numpy as np +import theano.tensor as tt +import theano + +from pymc3.distributions.distribution import ( + Distribution, Discrete, Continuous, draw_values, generate_samples) +from pymc3.distributions import transforms +from pymc3.distributions.dist_math import bound + +__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) + + 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(_Bounded, self).__init__( + shape=self._wrapped.shape, + dtype=self._wrapped.dtype, + testval=self._wrapped.testval, + defaults=defaults, + transform=self._wrapped.transform) + + def logp(self, value): + 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.') + samples = np.zeros(size, dtype=self.dtype).flatten() + i, n = 0, len(samples) + while i < len(samples): + sample = self._wrapped.random(point=point, size=n) + select = sample[np.logical_and(sample >= lower, sample <= upper)] + samples[i:(i + len(select))] = select[:] + i += len(select) + n -= len(select) + if size is not None: + return np.reshape(samples, size) + else: + return samples + + def random(self, point=None, size=None, repeat=None): + 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) + return generate_samples(self._random, lower, upper, point, + dist_shape=self.shape, + size=size) + elif self.lower is not None: + lower = draw_values([self.lower], point=point) + return generate_samples(self._random, lower, np.inf, point, + dist_shape=self.shape, + size=size) + else: + upper = draw_values([self.upper], point=point) + return generate_samples(self._random, -np.inf, upper, point, + dist_shape=self.shape, + size=size) + + +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 + + super(_DiscreteBounded, self).__init__( + distribution=distribution, lower=lower, upper=upper, + default=default, *args, **kwargs) + + +class _ContinuousBounded(_Bounded, Continuous): + R""" + An upper, lower or upper+lower bounded distribution + + 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 __init__(self, distribution, lower, upper, + transform='infer', *args, **kwargs): + dtype = kwargs.get('dtype', theano.config.floatX) + + if lower is not None: + lower = tt.as_tensor_variable(lower).astype(dtype) + if upper is not None: + upper = tt.as_tensor_variable(upper).astype(dtype) + + 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 + + super(_ContinuousBounded, self).__init__( + distribution=distribution, lower=lower, upper=upper, + transform=transform, default=default, *args, **kwargs) + + +class Bound(object): + R""" + Create a new upper, lower or upper+lower bounded distribution. + + The resulting distribution is not normalized anymore. This + is usually fine if the bounds are constants. If you need + truncated distributions, use `Bound` in combination with + a `pm.Potential` with the cumulative probability function. + + The bounds are inclusive for discrete distributions. + + Parameters + ---------- + distribution : pymc3 distribution + Distribution to be transformed into a bounded distribution. + lower : float or array like, optional + Lower bound of the distribution. + upper : float or array like, optional + Upper bound of the distribution. + + Example + ------- + # Bounded distribution can be defined before the model context + PositiveNormal = pm.Bound(pm.Normal, lower=0.0) + with pm.Model(): + par1 = PositiveNormal('par1', mu=0.0, sd=1.0, testval=1.0) + # or within the model context + NegativeNormal = pm.Bound(pm.Normal, upper=0.0) + par2 = NegativeNormal('par2', mu=0.0, sd=1.0, testval=1.0) + + # or you can define it implicitly within the model context + par3 = pm.Bound(pm.Normal, lower=-1.0, upper=1.0)( + 'par3', mu=0.0, sd=1.0, testval=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 + + if not issubclass(distribution, Distribution): + raise ValueError('"distribution" must be a Distribution subclass.') + + self.distribution = distribution + self.lower = lower + self.upper = upper + + def __call__(self, *args, **kwargs): + if 'observed' in kwargs: + raise ValueError('Observed Bound distributions are not supported. ' + 'If you want to model truncated data ' + 'you can use a pm.Potential in combination ' + 'with the cumulative probability function. See ' + 'pymc3/examples/censored_data.py for an example.') + first, args = args[0], args[1:] + + if issubclass(self.distribution, Continuous): + return _ContinuousBounded(first, self.distribution, + self.lower, self.upper, *args, **kwargs) + elif issubclass(self.distribution, Discrete): + return _DiscreteBounded(first, self.distribution, + self.lower, self.upper, *args, **kwargs) + else: + raise ValueError('Distribution is neither continuous nor discrete.') + + def dist(self, *args, **kwargs): + if issubclass(self.distribution, Continuous): + return _ContinuousBounded.dist( + self.distribution, self.lower, self.upper, *args, **kwargs) + + 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.') diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 94593b9c07..19186e16a8 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -21,7 +21,8 @@ bound, logpow, gammaln, betaln, std_cdf, i0, i1, alltrue_elemwise, SplineWrapper ) -from .distribution import Continuous, draw_values, generate_samples, Bound +from .distribution import Continuous, draw_values, generate_samples +from .bound import Bound __all__ = ['Uniform', 'Flat', 'Normal', 'Beta', 'Exponential', 'Laplace', 'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull', diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index d926d27e25..fdb4787118 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -6,12 +6,9 @@ from ..memoize import memoize from ..model import Model, get_named_nodes, FreeRV, ObservedRV from ..vartypes import string_types -from .dist_math import bound -# To avoid circular import for transform below -import pymc3 as pm -__all__ = ['DensityDist', 'Distribution', 'Continuous', 'Bound', - 'Discrete', 'NoDistribution', 'TensorType', 'draw_values'] +__all__ = ['DensityDist', 'Distribution', 'Continuous', 'Discrete', + 'NoDistribution', 'TensorType', 'draw_values'] class _Unpickling(object): @@ -396,126 +393,3 @@ def generate_samples(generator, *args, **kwargs): prefix_shape, *args, **kwargs) return reshape_sampled(samples, size, dist_shape) - - -class Bounded(Distribution): - R""" - An upper, lower or upper+lower bounded distribution - - 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 __init__(self, distribution, lower, upper, transform='infer', *args, **kwargs): - self.dist = distribution.dist(*args, **kwargs) - - self.__dict__.update(self.dist.__dict__) - self.__dict__.update(locals()) - - if hasattr(self.dist, 'mode'): - self.mode = self.dist.mode - - if transform == 'infer': - - default = self.dist.default() - - if not np.isinf(lower) and not np.isinf(upper): - self.transform = pm.distributions.transforms.interval(lower, upper) - if default <= lower or default >= upper: - self.testval = 0.5 * (upper + lower) - - if not np.isinf(lower) and np.isinf(upper): - self.transform = pm.distributions.transforms.lowerbound(lower) - if default <= lower: - self.testval = lower + 1 - - if np.isinf(lower) and not np.isinf(upper): - self.transform = pm.distributions.transforms.upperbound(upper) - if default >= upper: - self.testval = upper - 1 - - if issubclass(distribution, Discrete): - self.transform = None - - def _random(self, lower, upper, point=None, size=None): - samples = np.zeros(size).flatten() - i, n = 0, len(samples) - while i < len(samples): - sample = self.dist.random(point=point, size=n) - select = sample[np.logical_and(sample > lower, sample <= upper)] - samples[i:(i + len(select))] = select[:] - i += len(select) - n -= len(select) - if size is not None: - return np.reshape(samples, size) - else: - return samples - - def random(self, point=None, size=None, repeat=None): - lower, upper = draw_values([self.lower, self.upper], point=point) - return generate_samples(self._random, lower, upper, point, - dist_shape=self.shape, - size=size) - - def logp(self, value): - return bound(self.dist.logp(value), - value >= self.lower, value <= self.upper) - - -class Bound(object): - R""" - Creates a new upper, lower or upper+lower bounded distribution - - Parameters - ---------- - distribution : pymc3 distribution - Distribution to be transformed into a bounded distribution - lower : float (optional) - Lower bound of the distribution - upper : float (optional) - - Example - ------- - # Bounded distribution can be defined before the model context - PositiveNormal = pm.Bound(pm.Normal, lower=0.0) - with pm.Model(): - par1 = PositiveNormal('par1', mu=0.0, sd=1.0, testval=1.0) - # or within the model context - NegativeNormal = pm.Bound(pm.Normal, upper=0.0) - par2 = NegativeNormal('par2', mu=0.0, sd=1.0, testval=1.0) - - # or you can define it implicitly within the model context - par3 = pm.Bound(pm.Normal, lower=-1.0, upper=1.0)( - 'par3', mu=0.0, sd=1.0, testval=1.0) - """ - - def __init__(self, distribution, lower=-np.inf, upper=np.inf): - self.distribution = distribution - self.lower = lower - self.upper = upper - - def __call__(self, *args, **kwargs): - if 'observed' in kwargs: - raise ValueError('Observed Bound distributions are not allowed. ' - 'If you want to model truncated data ' - 'you can use a pm.Potential in combination ' - 'with the cumulative probability function. See ' - 'pymc3/examples/censored_data.py for an example.') - first, args = args[0], args[1:] - - return Bounded(first, self.distribution, self.lower, self.upper, - *args, **kwargs) - - def dist(self, *args, **kwargs): - return Bounded.dist(self.distribution, self.lower, self.upper, - *args, **kwargs) diff --git a/pymc3/distributions/transforms.py b/pymc3/distributions/transforms.py index 30f8e3a3b8..c813c44235 100644 --- a/pymc3/distributions/transforms.py +++ b/pymc3/distributions/transforms.py @@ -8,7 +8,7 @@ import numpy as np __all__ = ['transform', 'stick_breaking', 'logodds', 'interval', - 'lowerbound', 'upperbound', 'log', 'sum_to_1', 't_stick_breaking'] + 'lowerbound', 'upperbound', 'log', 'sum_to_1', 't_stick_breaking'] class Transform(object): @@ -33,6 +33,7 @@ def jacobian_det(self, x): raise NotImplementedError def apply(self, dist): + # avoid circular import return TransformedDistribution.dist(dist, self) def __str__(self): @@ -90,7 +91,7 @@ def backward(self, x): def forward(self, x): return tt.log(x) - + def forward_val(self, x, point=None): return self.forward(x) @@ -111,7 +112,7 @@ def backward(self, x): def forward(self, x): return logit(x) - + def forward_val(self, x, point=None): return self.forward(x) @@ -324,6 +325,6 @@ def forward(self, y): def forward_val(self, x, point=None): return self.forward(x) - + def jacobian_det(self, y): return tt.sum(y[self.diag_idxs]) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index ca012495ce..0edccc3e49 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -18,7 +18,7 @@ from ..distributions import continuous from pymc3.theanof import floatX from numpy import array, inf, log, exp -from numpy.testing import assert_almost_equal +from numpy.testing import assert_almost_equal, assert_equal import numpy.random as nr import numpy as np import pytest @@ -569,7 +569,8 @@ def test_binomial(self): self.pymc3_matches_scipy(Binomial, Nat, {'n': NatSmall, 'p': Unit}, lambda value, n, p: sp.binom.logpmf(value, n, p)) - @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") # Too lazy to propagate decimal parameter through the whole chain of deps + # Too lazy to propagate decimal parameter through the whole chain of deps + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_beta_binomial(self): self.checkd(BetaBinomial, Nat, {'alpha': Rplus, 'beta': Rplus, 'n': NatSmall}) @@ -597,16 +598,19 @@ def test_constantdist(self): self.pymc3_matches_scipy(Constant, I, {'c': I}, lambda value, c: np.log(c == value)) - @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") # Too lazy to propagate decimal parameter through the whole chain of deps + # Too lazy to propagate decimal parameter through the whole chain of deps + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_zeroinflatedpoisson(self): self.checkd(ZeroInflatedPoisson, Nat, {'theta': Rplus, 'psi': Unit}) - @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") # Too lazy to propagate decimal parameter through the whole chain of deps + # Too lazy to propagate decimal parameter through the whole chain of deps + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_zeroinflatednegativebinomial(self): self.checkd(ZeroInflatedNegativeBinomial, Nat, {'mu': Rplusbig, 'alpha': Rplusbig, 'psi': Unit}) - @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") # Too lazy to propagate decimal parameter through the whole chain of deps + # Too lazy to propagate decimal parameter through the whole chain of deps + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_zeroinflatedbinomial(self): self.checkd(ZeroInflatedBinomial, Nat, {'n': NatSmall, 'p': Unit, 'psi': Unit}) @@ -850,6 +854,68 @@ def ref_pdf(value): self.pymc3_matches_scipy(TestedInterpolated, R, {}, ref_pdf) +def test_bound(): + np.random.seed(42) + UnboundNormal = Bound(Normal) + dist = UnboundNormal.dist(mu=0, sd=1) + assert dist.transform is None + assert dist.default() == 0. + assert isinstance(dist.random(), np.ndarray) + + LowerNormal = Bound(Normal, lower=1) + dist = LowerNormal.dist(mu=0, sd=1) + assert dist.logp(0).eval() == -np.inf + assert dist.default() > 1 + assert dist.transform is not None + assert np.all(dist.random() > 1) + + UpperNormal = Bound(Normal, upper=-1) + dist = UpperNormal.dist(mu=0, sd=1) + assert dist.logp(-0.5).eval() == -np.inf + assert dist.default() < -1 + 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, sd=1, shape=2) + assert_equal(dist.logp([0.5, 3.5]).eval(), -np.array([np.inf, np.inf])) + assert_equal(dist.default(), np.array([1.5, 2.5])) + 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', shape=2) + assert_equal(a.tag.test_value, np.array([1.5, 2.5])) + + lower = tt.vector('lower') + lower.tag.test_value = np.array([1, 2]).astype(theano.config.floatX) + upper = 3 + ArrayNormal = Bound(Normal, lower=lower, upper=upper) + dist = ArrayNormal.dist(mu=0, sd=1, shape=2) + logp = dist.logp([0.5, 3.5]).eval({lower: lower.tag.test_value}) + assert_equal(logp, -np.array([np.inf, np.inf])) + assert_equal(dist.default(), np.array([2, 2.5])) + assert dist.transform is not None + + with Model(): + a = ArrayNormal('c', shape=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 + + def test_repr_latex_(): with Model(): x0 = Binomial('Discrete', p=.5, n=10) diff --git a/pymc3/tests/test_examples.py b/pymc3/tests/test_examples.py index 415d96a427..79d3736c79 100644 --- a/pymc3/tests/test_examples.py +++ b/pymc3/tests/test_examples.py @@ -6,6 +6,7 @@ import theano.tensor as tt import pytest import theano +from pymc3.theanof import floatX from .helpers import SeededTest @@ -43,20 +44,14 @@ def build_model(self): with pm.Model() as model: effects = pm.Normal('effects', mu=0, tau=100. ** -2, shape=len(P.columns)) - p = tt.nnet.sigmoid(tt.dot(np.array(P), effects)) - pm.Bernoulli('s', p, observed=np.array(data.switch)) + p = tt.nnet.sigmoid(tt.dot(floatX(np.array(P)), effects)) + pm.Bernoulli('s', p, observed=floatX(np.array(data.switch))) return model def test_run(self): model = self.build_model() with model: - # move the chain to the MAP which should be a good starting point - start = pm.find_MAP() - H = model.fastd2logp() # find a good orientation using the hessian at the MAP - h = H(start) - - step = pm.HamiltonianMC(model.vars, h) - pm.sample(50, step=step, start=start) + pm.sample(50, tune=50, n_init=1000) class TestARM12_6(SeededTest):