Fix bug in Bound

Author: aseyboldt

pymc3/distributions/distribution.py

@@ -444,17 +444,20 @@ def __init__(self, distribution, lower, upper,
             else:
                 transform = transforms.lowerbound(lower)
                 default = lower + 1
+        else:
+            default = None
 
-        # We only change logp and testval for
-        # discrete distributions
+        # We don't use transformations for dicrete variables
         if issubclass(distribution, Discrete):
             transform = None
-            if default is not None:
-                default = default.astype(self.dist.default.type())
 
+        kwargs['transform'] = transform
         self._wrapped = distribution.dist(*args, **kwargs)
         self._default = default
 
+        if issubclass(distribution, Discrete) and default is not None:
+            default = default.astype(str(self._wrapped.default().dtype))
+
         if default is None:
             defaults = self._wrapped.defaults
             for name in defaults:
@@ -470,16 +473,16 @@ def __init__(self, distribution, lower, upper,
             transform=self._wrapped.transform)
 
     def _random(self, lower, upper, point=None, size=None):
-        if lower is None:
-            lower = -np.inf
-        if upper is None:
-            upper = np.inf
-
-        samples = np.zeros(size).flatten()
+        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)]
+            select = sample[np.logical_and(sample >= lower, sample <= upper)]
             samples[i:(i + len(select))] = select[:]
             i += len(select)
             n -= len(select)
@@ -489,18 +492,31 @@ def _random(self, lower, upper, point=None, size=None):
             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)
+        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)
 
     def logp(self, value):
         logp = self._wrapped.logp(value)
         bounds = []
         if self.lower is not None:
-            bounds.append(value > self.lower)
+            bounds.append(value >= self.lower)
         if self.upper is not None:
-            bounds.append(value < self.upper)
+            bounds.append(value <= self.upper)
         if len(bounds) > 0:
             return bound(logp, *bounds)
         else:
@@ -516,13 +532,16 @@ class Bound(object):
     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
+        Distribution to be transformed into a bounded distribution.
     lower : float or array like, optional
-        Lower bound of the distribution
+        Lower bound of the distribution.
     upper : float or array like, optional
+        Upper bound of the distribution.
 
     Example
     -------

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

pymc3/tests/test_distributions.py

@@ -855,26 +855,44 @@ def ref_pdf(value):
 
 
 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(-1).eval() == -np.inf
+    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)
-    assert_equal(dist.logp([1, 2]).eval(), -np.array([np.inf, np.inf]))
+    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]))
 
+    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
+
 
 def test_repr_latex_():
     with Model():