Skip to content

quick solution to Bounded sampling of arbitrary shapes #1069

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 3 commits into from
Closed

quick solution to Bounded sampling of arbitrary shapes #1069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b50c566
quick solution to Bounded sampling of arbitrary shapes
brandonwillard Apr 13, 2016
8fbfe36
added a iid option to the Bound sampler
brandonwillard Apr 20, 2016
a41fd67
changes test sizes for the bounded distribution cases
brandonwillard Apr 26, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 46 additions & 9 deletions pymc3/distributions/continuous.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -991,24 +991,61 @@ def __init__(self, distribution, lower, upper, *args, **kwargs):
if hasattr(self.dist, 'mode'):
self.mode = self.dist.mode

def _random(self, lower, upper, point=None, size=None):
samples = np.zeros(size).flatten()
def _random(self, lower, upper, point=None, iid=False, **kwargs):
r"""Rejection approach to truncated sampling.

Samples from `self.dist` (i.e. the underlying distribution)
and accepts samples, :math:`x`, that satisfy
.. math::

\mathtt{lower} < x \leq \mathtt{upper}.

Parameters
----------
lower : numpy.ndarray
Lower bound. Must be comparable to `self.dist.shape`.
upper : numpy.ndarray
Upper bound. Must be comparable to `self.dist.shape`.
point : dict or None
Conditional parameter values?
iid : bool
Are the samples independent and identically distributed?

Returns
-------
numpy.ndarray
A single truncated sample from `self.dist`.

Notes
-----
Assumes that `self.dist` is a `self.shape`-tensor consists of
i.i.d. random variables.
"""
samples = np.zeros(shape=self.shape).flatten()
shape_len = self.dist.shape.prod()
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[:]
q, r = divmod(n, shape_len)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just realized that this was left over from an attempt to make _random behave as it used to; specifically, when it took the size parameter. My initial implementation had samples take the shape size + tuple(self.shape), then code would produce the desired result from this method alone. However, after looking into the generating_samples function, I noticed that it takes steps to produce replications with shape given by size.

sample_size = q + min(1, r)
sample = self.dist.random(point=point, size=(sample_size,))
accepted_ind = np.logical_and(sample > lower, sample <= upper)

if not iid and not np.all(accepted_ind):
continue

select = sample[accepted_ind]
sel_len = min(n, len(select))
samples[i:(i+sel_len)] = select[0:sel_len]
i += len(select)
n -= len(select)
if size is not None:
return np.reshape(samples, size)
else:
return samples

return np.reshape(samples, self.shape)

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,
broadcast_shape=tuple(self.shape),
size=size)

def logp(self, value):
Expand Down
57 changes: 40 additions & 17 deletions pymc3/tests/test_distributions_random.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,9 +94,9 @@ def check_dist(dist_case, test_cases, shape=None):
dist, dist_kwargs = dist_case
with Model():
if shape is None:
rv = dist(dist.__name__, transform=None, **dist_kwargs)
rv = dist(getattr(dist, "__name__", "none"), transform=None, **dist_kwargs)
else:
rv = dist(dist.__name__, shape=shape, transform=None,
rv = dist(getattr(dist, "__name__", "none"), shape=shape, transform=None,
**dist_kwargs)
for size, expected in test_cases:
check_shape(rv, size=size, expected=expected)
Expand All @@ -119,13 +119,18 @@ def check_shape(rv, size=None, expected=None):
@attr('scalar_parameter_shape')
class ScalarParameterShape(unittest.TestCase):

def check(self, dist, **kwargs):
test_cases = [(None, (1,)), (5, (5,)), ((4, 5), (4, 5))]
def check(self, dist, test_cases=None, **kwargs):
if test_cases is None:
test_cases = [(None, (1,)), (5, (5,)), ((4, 5), (4, 5))]
check_dist((dist, kwargs), test_cases)

def test_normal(self):
self.check(Normal, mu=0., tau=1.)

def test_bounded(self):
BoundedNormal = Bound(Normal, upper=0)
self.check(BoundedNormal, mu=0., tau=1.)

def test_uniform(self):
self.check(Uniform, lower=0., upper=1.)

Expand Down Expand Up @@ -214,14 +219,18 @@ def test_categorical(self):
@attr('scalar_shape')
class ScalarShape(unittest.TestCase):

def check(self, dist, **kwargs):
n = 10
test_cases = [(None, (n,)), (5, (5, n,)), ((4, 5), (4, 5, n,))]
check_dist((dist, kwargs), test_cases, n)
def check(self, dist, _n=10, test_cases=None, **kwargs):
if test_cases is None:
test_cases = [(None, (_n,)), (5, (5, _n,)), ((4, 5), (4, 5, _n,))]
check_dist((dist, kwargs), test_cases, _n)

def test_normal(self):
self.check(Normal, mu=0., tau=1.)

def test_bounded(self):
BoundedNormal = Bound(Normal, upper=0)
self.check(BoundedNormal, _n=3, mu=0., tau=1.)

def test_uniform(self):
self.check(Uniform, lower=0., upper=1.)

Expand Down Expand Up @@ -314,14 +323,21 @@ def setUp(self):
self.zeros = np.zeros(self.n)
self.ones = np.ones(self.n)

def check(self, dist, **kwargs):
n = self.n
test_cases = [(None, (n,)), (5, (5, n,)), ((4, 5), (4, 5, n,))]
check_dist((dist, kwargs), test_cases, n)
def check(self, dist, _n=None, test_cases=None, **kwargs):
if _n is None:
_n = self.n
if test_cases is None:
test_cases = [(None, (_n,)), (5, (5, _n,)), ((4, 5), (4, 5, _n,))]
check_dist((dist, kwargs), test_cases, _n)

def test_normal(self):
self.check(Normal, mu=self.zeros, tau=self.ones)

def test_bounded(self):
BoundedNormal = Bound(Normal, upper=0)
n = 2
self.check(BoundedNormal, _n=n, mu=np.zeros(n), tau=np.ones(n))

def test_uniform(self):
self.check(Uniform, lower=self.zeros, upper=self.ones)

Expand Down Expand Up @@ -425,16 +441,23 @@ def setUp(self):
self.zeros = np.zeros(self.n)
self.ones = np.ones(self.n)

def check(self, dist, **kwargs):
n = self.n
shape = (2*n, n)
test_cases = [(None, shape), (5, (5,) + shape),
((4, 5), (4, 5) + shape)]
def check(self, dist, _n=None, test_cases=None, **kwargs):
if _n is None:
_n = self.n
shape = (2*_n, _n)
if test_cases is None:
test_cases = [(None, shape), (5, (5,) + shape),
((4, 5), (4, 5) + shape)]
check_dist((dist, kwargs), test_cases, shape)

def test_normal(self):
self.check(Normal, mu=self.zeros, tau=self.ones)

def test_bounded(self):
BoundedNormal = Bound(Normal, upper=0)
n = 2
self.check(BoundedNormal, _n=n, mu=np.zeros(n), tau=np.ones(n))

def test_uniform(self):
self.check(Uniform, lower=self.zeros, upper=self.ones)

Expand Down