From 3a9a2955b975569d35458b37e126741fe1ceee39 Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Wed, 14 Nov 2018 23:06:46 +0100 Subject: [PATCH 1/4] Fix for #3225. Made Triangular `c` attribute be handled consistently with scipy.stats. Added test and updated example code. --- RELEASE-NOTES.md | 1 + pymc3/distributions/continuous.py | 18 +++++++++++------- pymc3/tests/test_distributions_random.py | 9 ++++++++- 3 files changed, 20 insertions(+), 8 deletions(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 31c5bcb175..a493e3135a 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -13,6 +13,7 @@ ### Maintenance +- Fixed Triangular distribution `c` attribute handling in `random` and updated sample codes for consistency (#3225) - Renamed `sample_ppc()` and `sample_ppc_w()` to `sample_posterior_predictive()` and `sample_posterior_predictive_w()`, respectively. - Refactor SMC and properly compute marginal likelihood (#3124) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 984ca3e577..f0034e51d4 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -3262,13 +3262,15 @@ class Triangular(BoundedContinuous): plt.style.use('seaborn-darkgrid') x = np.linspace(-2, 10, 500) lowers = [0., -1, 2] - cs = [0.5, 0.5, 0.75] - uppers = [4., 2, 6] - for lower, c_, upper_ in zip(lowers, cs, uppers): - pdf = st.triang.pdf(x, loc=lower, c=c_, scale=upper_) + cs = [2., 0., 6.5] + uppers = [4., 1, 8] + for lower, c, upper in zip(lowers, cs, uppers): + scale = upper - lower + c_ = (c - lower) / scale + pdf = st.triang.pdf(x, loc=lower, c=c_, scale=scale) plt.plot(x, pdf, label='lower = {}, c = {}, upper = {}'.format(lower, - lower + upper_ * c_, - lower + upper_)) + c, + upper)) plt.xlabel('x', fontsize=12) plt.ylabel('f(x)', fontsize=12) plt.legend(loc=1) @@ -3318,7 +3320,9 @@ def random(self, point=None, size=None): """ c, lower, upper = draw_values([self.c, self.lower, self.upper], point=point, size=size) - return generate_samples(stats.triang.rvs, c=c-lower, loc=lower, scale=upper-lower, + scale = upper - lower + c_ = (c - lower) / scale + return generate_samples(stats.triang.rvs, c=c_, loc=lower, scale=scale, size=size, dist_shape=self.shape, random_state=None) def logp(self, value): diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 9bc2841624..0e354536c6 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -13,7 +13,7 @@ from pymc3.distributions.distribution import draw_values from .helpers import SeededTest from .test_distributions import ( - build_model, Domain, product, R, Rplus, Rplusbig, Rplusdunif, + build_model, Domain, product, R, Rplus, Rplusbig, Runif, Rplusdunif, Unit, Nat, NatSmall, I, Simplex, Vector, PdMatrix, PdMatrixChol, PdMatrixCholUpper, RealMatrix, RandomPdMatrix ) @@ -518,6 +518,13 @@ def ref_rand(size, mu, kappa): return st.vonmises.rvs(size=size, loc=mu, kappa=kappa) pymc3_random(pm.VonMises, {'mu': R, 'kappa': Rplus}, ref_rand=ref_rand) + def test_triangular(self): + def ref_rand(size, lower, upper, c): + scale = upper - lower + c_ = (c - lower) / scale + return st.triang.rvs(size=size, loc=lower, scale=scale, c=c_) + pymc3_random(pm.Triangular, {'lower': Runif, 'upper': Runif + 3, 'c': Runif + 1}, ref_rand=ref_rand) + def test_flat(self): with pm.Model(): f = pm.Flat('f') From 6cefd1751bec140e7337aae644d65189558ff05a Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Tue, 27 Nov 2018 17:15:18 +0100 Subject: [PATCH 2/4] Fix for #3210 which uses a completely different approach than PR #3214. It uses a context manager inside `draw_values` that makes all the values drawn from `TensorVariables` or `MultiObservedRV`s available to nested calls of the original call to `draw_values`. It is partly inspired by how Edward2 approaches the problem of forward sampling. Ed2 tensors fix a `_values` attribute after they first call `sample` and then only return that. They can do it because of their functional scheme, where the entire graph is recreated each time the generative function is called. Our object oriented paradigm cannot set a fixed _values, it has to know it is in the context of a single `draw_values` call. That is why I opted for context managers to store the drawn values. --- pymc3/distributions/distribution.py | 311 ++++++++++++++------- pymc3/model.py | 10 + pymc3/tests/test_random.py | 58 ++++ pymc3/tests/test_step.py | 400 ++++++++++++++-------------- 4 files changed, 482 insertions(+), 297 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 5783609765..e5681ea34c 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -1,4 +1,4 @@ -import collections +import six import numbers import numpy as np @@ -8,7 +8,7 @@ from ..memoize import memoize from ..model import ( Model, get_named_nodes_and_relations, FreeRV, - ObservedRV, MultiObservedRV + ObservedRV, MultiObservedRV, Context, InitContextMeta ) from ..vartypes import string_types @@ -214,6 +214,48 @@ def random(self, *args, **kwargs): "Define a custom random method and pass it as kwarg random") +class _DrawValuesContext(six.with_metaclass(InitContextMeta, Context)): + """ A context manager class used while drawing values with draw_values + """ + + def __new__(cls, *args, **kwargs): + # resolves the parent instance + instance = super(_DrawValuesContext, cls).__new__(cls) + if cls.get_contexts(): + potencial_parent = cls.get_contexts()[-1] + # We have to make sure that the context is a _DrawValuesContext + # and not a Model + if isinstance(potencial_parent, cls): + instance._parent = potencial_parent + else: + instance._parent = None + else: + instance._parent = None + return instance + + def __init__(self): + if self.parent is not None: + # All _DrawValuesContext instances that are in the context of + # another _DrawValuesContext will share the reference to the + # drawn_vars dictionary. This means that separate branches + # in the nested _DrawValuesContext context tree will see the + # same drawn values + self.drawn_vars = self.parent.drawn_vars + else: + self.drawn_vars = dict() + + @property + def parent(self): + return self._parent + + +def is_fast_drawable(var): + return isinstance(var, (numbers.Number, + np.ndarray, + tt.TensorConstant, + tt.sharedvar.SharedVariable)) + + def draw_values(params, point=None, size=None): """ Draw (fix) parameter values. Handles a number of cases: @@ -232,97 +274,134 @@ def draw_values(params, point=None, size=None): b) are *RVs with a random method """ - # Distribution parameters may be nodes which have named node-inputs - # specified in the point. Need to find the node-inputs, their - # parents and children to replace them. - leaf_nodes = {} - named_nodes_parents = {} - named_nodes_children = {} - for param in params: - if hasattr(param, 'name'): - # Get the named nodes under the `param` node - nn, nnp, nnc = get_named_nodes_and_relations(param) - leaf_nodes.update(nn) - # Update the discovered parental relationships - for k in nnp.keys(): - if k not in named_nodes_parents.keys(): - named_nodes_parents[k] = nnp[k] - else: - named_nodes_parents[k].update(nnp[k]) - # Update the discovered child relationships - for k in nnc.keys(): - if k not in named_nodes_children.keys(): - named_nodes_children[k] = nnc[k] - else: - named_nodes_children[k].update(nnc[k]) - - # Init givens and the stack of nodes to try to `_draw_value` from - givens = {} - stored = set() # Some nodes - stack = list(leaf_nodes.values()) # A queue would be more appropriate - while stack: - next_ = stack.pop(0) - if next_ in stored: - # If the node already has a givens value, skip it - continue - elif isinstance(next_, (tt.TensorConstant, - tt.sharedvar.SharedVariable)): - # If the node is a theano.tensor.TensorConstant or a - # theano.tensor.sharedvar.SharedVariable, its value will be - # available automatically in _compile_theano_function so - # we can skip it. Furthermore, if this node was treated as a - # TensorVariable that should be compiled by theano in - # _compile_theano_function, it would raise a `TypeError: - # ('Constants not allowed in param list', ...)` for - # TensorConstant, and a `TypeError: Cannot use a shared - # variable (...) as explicit input` for SharedVariable. - stored.add(next_.name) - continue - else: - # If the node does not have a givens value, try to draw it. - # The named node's children givens values must also be taken - # into account. - children = named_nodes_children[next_] - temp_givens = [givens[k] for k in givens if k in children] - try: - # This may fail for autotransformed RVs, which don't - # have the random method - givens[next_.name] = (next_, _draw_value(next_, - point=point, - givens=temp_givens, - size=size)) - stored.add(next_.name) - except theano.gof.fg.MissingInputError: - # The node failed, so we must add the node's parents to - # the stack of nodes to try to draw from. We exclude the - # nodes in the `params` list. - stack.extend([node for node in named_nodes_parents[next_] - if node is not None and - node.name not in stored and - node not in params]) - - # the below makes sure the graph is evaluated in order - # test_distributions_random::TestDrawValues::test_draw_order fails without it - params = dict(enumerate(params)) # some nodes are not hashable - evaluated = {} - to_eval = set() - missing_inputs = set(params) - while to_eval or missing_inputs: - if to_eval == missing_inputs: - raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval])) - to_eval = set(missing_inputs) - missing_inputs = set() - for param_idx in to_eval: - param = params[param_idx] - if hasattr(param, 'name') and param.name in givens: - evaluated[param_idx] = givens[param.name][1] + # Get fast drawable values (i.e. things in point or numbers, arrays, + # constants or shares, or things that were already drawn in related + # contexts) + if point is None: + point = {} + with _DrawValuesContext() as context: + params = dict(enumerate(params)) + drawn = context.drawn_vars + evaluated = {} + symbolic_params = [] + for i, p in params.items(): + # If the param is fast drawable, then draw the value immediately + if is_fast_drawable(p): + v = _draw_value(p, point=point, size=size) + evaluated[i] = v + continue + + name = getattr(p, 'name', None) + if p in drawn: + # param was drawn in related contexts + v = drawn[p] + evaluated[i] = v + elif name is not None and name in point: + # param.name is in point + v = point[name] + evaluated[i] = drawn[p] = v else: - try: # might evaluate in a bad order, - evaluated[param_idx] = _draw_value(param, point=point, givens=givens.values(), size=size) - if isinstance(param, collections.Hashable) and named_nodes_parents.get(param): - givens[param.name] = (param, evaluated[param_idx]) + # param still needs to be drawn + symbolic_params.append((i, p)) + + if not symbolic_params: + # We only need to enforce the correct order if there are symbolic + # params that could be drawn in variable order + return [evaluated[i] for i in params] + + # Distribution parameters may be nodes which have named node-inputs + # specified in the point. Need to find the node-inputs, their + # parents and children to replace them. + leaf_nodes = {} + named_nodes_parents = {} + named_nodes_children = {} + for _, param in symbolic_params: + if hasattr(param, 'name'): + # Get the named nodes under the `param` node + nn, nnp, nnc = get_named_nodes_and_relations(param) + leaf_nodes.update(nn) + # Update the discovered parental relationships + for k in nnp.keys(): + if k not in named_nodes_parents.keys(): + named_nodes_parents[k] = nnp[k] + else: + named_nodes_parents[k].update(nnp[k]) + # Update the discovered child relationships + for k in nnc.keys(): + if k not in named_nodes_children.keys(): + named_nodes_children[k] = nnc[k] + else: + named_nodes_children[k].update(nnc[k]) + + # Init givens and the stack of nodes to try to `_draw_value` from + givens = {p.name: (p, v) for p, v in drawn.items() + if getattr(p, 'name', None) is not None} + stack = list(leaf_nodes.values()) # A queue would be more appropriate + while stack: + next_ = stack.pop(0) + if next_ in drawn: + # If the node already has a givens value, skip it + continue + elif isinstance(next_, (tt.TensorConstant, + tt.sharedvar.SharedVariable)): + # If the node is a theano.tensor.TensorConstant or a + # theano.tensor.sharedvar.SharedVariable, its value will be + # available automatically in _compile_theano_function so + # we can skip it. Furthermore, if this node was treated as a + # TensorVariable that should be compiled by theano in + # _compile_theano_function, it would raise a `TypeError: + # ('Constants not allowed in param list', ...)` for + # TensorConstant, and a `TypeError: Cannot use a shared + # variable (...) as explicit input` for SharedVariable. + continue + else: + # If the node does not have a givens value, try to draw it. + # The named node's children givens values must also be taken + # into account. + children = named_nodes_children[next_] + temp_givens = [givens[k] for k in givens if k in children] + try: + # This may fail for autotransformed RVs, which don't + # have the random method + value = _draw_value(next_, + point=point, + givens=temp_givens, + size=size) + givens[next_.name] = (next_, value) + drawn[next_] = value except theano.gof.fg.MissingInputError: - missing_inputs.add(param_idx) + # The node failed, so we must add the node's parents to + # the stack of nodes to try to draw from. We exclude the + # nodes in the `params` list. + stack.extend([node for node in named_nodes_parents[next_] + if node is not None and + node.name not in drawn and + node not in params]) + + # the below makes sure the graph is evaluated in order + # test_distributions_random::TestDrawValues::test_draw_order fails without it + # The remaining params that must be drawn are all hashable + to_eval = set() + missing_inputs = set([j for j, p in symbolic_params]) + while to_eval or missing_inputs: + if to_eval == missing_inputs: + raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval])) + to_eval = set(missing_inputs) + missing_inputs = set() + for param_idx in to_eval: + param = params[param_idx] + if param in drawn: + evaluated[param_idx] = drawn[param] + else: + try: # might evaluate in a bad order, + value = _draw_value(param, + point=point, + givens=givens.values(), + size=size) + evaluated[param_idx] = drawn[param] = value + givens[param.name] = (param, value) + except theano.gof.fg.MissingInputError: + missing_inputs.add(param_idx) return [evaluated[j] for j in params] # set the order back @@ -400,8 +479,16 @@ def _draw_value(param, point=None, givens=None, size=None): # reset shape to account for shape changes # with theano.shared inputs dist_tmp.shape = np.array([]) - val = dist_tmp.random(point=point, size=None) - dist_tmp.shape = val.shape + val = np.atleast_1d(dist_tmp.random(point=point, + size=None)) + # Sometimes point may change the size of val but not the + # distribution's shape + if point and size is not None: + temp_size = np.atleast_1d(size) + if all(val.shape[:len(temp_size)] == temp_size): + dist_tmp.shape = val.shape[len(temp_size):] + else: + dist_tmp.shape = val.shape return dist_tmp.random(point=point, size=size) else: return param.distribution.random(point=point, size=size) @@ -411,10 +498,24 @@ def _draw_value(param, point=None, givens=None, size=None): else: variables = values = [] func = _compile_theano_function(param, variables) - if size and values and not all(var.dshape == val.shape for var, val in zip(variables, values)): - return np.array([func(*v) for v in zip(*values)]) + if size is not None: + size = np.atleast_1d(size) + dshaped_variables = all((hasattr(var, 'dshape') + for var in variables)) + if (values and dshaped_variables and + not all(var.dshape == getattr(val, 'shape', tuple()) + for var, val in zip(variables, values))): + output = np.array([func(*v) for v in zip(*values)]) + elif (size is not None and any((val.ndim > var.ndim) + for var, val in zip(variables, values))): + output = np.array([func(*v) for v in zip(*values)]) else: - return func(*values) + output = func(*values) + return output + print(param, + type(param), + isinstance(param, tt.TensorVariable), + isinstance(param, (tt.TensorVariable, MultiObservedRV))) raise ValueError('Unexpected type in draw_value: %s' % type(param)) @@ -499,6 +600,20 @@ def generate_samples(generator, *args, **kwargs): samples = generator(size=broadcast_shape, *args, **kwargs) elif dist_shape == broadcast_shape: samples = generator(size=size_tup + dist_shape, *args, **kwargs) + elif len(dist_shape) == 0 and size_tup and broadcast_shape[:len(size_tup)] == size_tup: + # Input's dist_shape is scalar, but it has size repetitions. + # So now the size matches but we have to manually broadcast to + # the right dist_shape + samples = [generator(*args, **kwargs)] + if samples[0].shape == broadcast_shape: + samples = samples[0] + else: + suffix = broadcast_shape[len(size_tup):] + dist_shape + samples.extend([generator(*args, **kwargs). + reshape(broadcast_shape)[..., np.newaxis] + for _ in range(np.prod(suffix, + dtype=int) - 1)]) + samples = np.hstack(samples).reshape(size_tup + suffix) else: samples = None # Args have been broadcast correctly, can just ask for the right shape out @@ -515,9 +630,11 @@ def generate_samples(generator, *args, **kwargs): if samples is None: raise TypeError('''Attempted to generate values with incompatible shapes: size: {size} + size_tup: {size_tup} + broadcast_shape[:len(size_tup)] == size_tup: {test} dist_shape: {dist_shape} broadcast_shape: {broadcast_shape} - '''.format(size=size, dist_shape=dist_shape, broadcast_shape=broadcast_shape)) + '''.format(size=size, size_tup=size_tup, dist_shape=dist_shape, broadcast_shape=broadcast_shape, test=broadcast_shape[:len(size_tup)] == size_tup)) # reshape samples here if samples.shape[0] == 1 and size == 1: diff --git a/pymc3/model.py b/pymc3/model.py index 9d705de25a..c9ec85fac6 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -1386,6 +1386,16 @@ def __init__(self, name, data, distribution, total_size=None, model=None): self.distribution = distribution self.scaling = _get_scaling(total_size, self.logp_elemwiset.shape, self.logp_elemwiset.ndim) + # Make hashable by id for draw_values + def __hash__(self): + return id(self) + + def __eq__(self, other): + return self.id == other.id + + def __ne__(self, other): + return not self == other + def _walk_up_rv(rv): """Walk up theano graph to get inputs for deterministic RV.""" diff --git a/pymc3/tests/test_random.py b/pymc3/tests/test_random.py index f61fe3b4a7..1eed3dcd9b 100644 --- a/pymc3/tests/test_random.py +++ b/pymc3/tests/test_random.py @@ -1,11 +1,13 @@ import pymc3 as pm import numpy as np +from numpy import random as nr import numpy.testing as npt import pytest import theano.tensor as tt import theano from pymc3.distributions.distribution import _draw_value, draw_values +from .helpers import SeededTest def test_draw_value(): @@ -88,3 +90,59 @@ def test_dep_vars(self): assert all([np.all(val1 != val2), np.all(val1 != val3), np.all(val1 != val4), np.all(val2 != val3), np.all(val2 != val4), np.all(val3 != val4)]) + + +class TestJointDistributionDrawValues(SeededTest): + def test_joint_distribution(self): + with pm.Model() as model: + a = pm.Normal('a', mu=0, sd=100) + b = pm.Normal('b', mu=a, sd=1e-8) + c = pm.Normal('c', mu=a, sd=1e-8) + d = pm.Deterministic('d', b + c) + + # Expected RVs + N = 1000 + norm = np.random.randn(3, N) + eA = norm[0] * 100 + eB = eA + norm[1] * 1e-8 + eC = eA + norm[2] * 1e-8 + eD = eB + eC + + # Drawn RVs + nr.seed(self.random_seed) +# A, B, C, D = list(zip(*[draw_values([a, b, c, d]) for i in range(N)])) + A, B, C, D = draw_values([a, b, c, d], size=N) + A = np.array(A).flatten() + B = np.array(B).flatten() + C = np.array(C).flatten() + D = np.array(D).flatten() + + # Assert that the drawn samples match the expected values + assert np.allclose(eA, A) + assert np.allclose(eB, B) + assert np.allclose(eC, C) + assert np.allclose(eD, D) + + # Assert that A, B and C have the expected difference + assert np.all(np.abs(A - B) < 1e-6) + assert np.all(np.abs(A - C) < 1e-6) + assert np.all(np.abs(B - C) < 1e-6) + + # Marginal draws + mA = np.array([draw_values([a]) for i in range(N)]).flatten() + mB = np.array([draw_values([b]) for i in range(N)]).flatten() + mC = np.array([draw_values([c]) for i in range(N)]).flatten() + # Also test the with model context of draw_values + with model: + mD = np.array([draw_values([d]) for i in range(N)]).flatten() + + # Assert that the marginal distributions have different sample values + assert not np.all(np.abs(B - mB) < 1e-2) + assert not np.all(np.abs(C - mC) < 1e-2) + assert not np.all(np.abs(D - mD) < 1e-2) + + # Assert that the marginal distributions do not have high cross + # correlation + assert np.abs(np.corrcoef(mA, mB)[0, 1]) < 0.1 + assert np.abs(np.corrcoef(mA, mC)[0, 1]) < 0.1 + assert np.abs(np.corrcoef(mB, mC)[0, 1]) < 0.1 diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py index 4b1cfe17b8..28b1814dcb 100644 --- a/pymc3/tests/test_step.py +++ b/pymc3/tests/test_step.py @@ -467,206 +467,206 @@ class TestStepMethods(object): # yield test doesn't work subclassing object ), SMC: np.array( [ - 0.02957155, - -0.10448987, - 1.83330321, - 0.45125892, - 0.08847866, - 1.42446262, - -0.50373428, - 1.7926401, - -0.43595123, - 1.29051484, - 1.48825784, - 0.58940796, - 0.39452458, - 0.6458593, - 1.41254985, - 0.33640003, - -0.82261885, - 0.19660577, - -1.31652677, - -0.43394062, - 0.47844794, - 0.32073809, - 0.43811662, - 0.15053655, - 0.36644253, - 0.67791414, - -0.02508511, - 0.7298907, - 1.98694909, - 1.7186239, - 0.36670138, - 1.31585297, - 1.0010273, - 1.03616273, - 0.14999443, - 0.27743457, - -0.01086402, - 1.49293981, - 1.56440459, - -0.05856325, - 1.21029567, - 0.3526936, - 0.72190671, - 0.06857266, - 2.5848243, - 0.03416341, - 0.11731779, - 0.71426943, - -0.01581274, - 0.31147419, - 0.59718523, - 0.25670663, - 0.71827237, - 0.19043905, - 1.32676133, - -0.02546965, - 1.56694401, - -0.27818479, - -0.32227988, - -0.58671167, - 0.18231422, - -0.51437682, - 0.60519212, - 2.01615164, - 1.25931896, - 1.0957265, - 2.35609054, - -1.03617366, - -0.38976571, - 1.2889846, - 0.39140397, - -0.88146889, - 0.11897506, - 0.68374636, - 0.35969604, - 1.55791486, - 0.60414685, - 0.76464213, - 0.5734933, - 0.55888177, - 1.14165163, - 0.29172266, - 1.0656623, - 1.49059082, - -0.78259634, - 1.39116132, - -1.30272835, - 1.47197121, - -0.44181174, - 0.61235063, - 0.06672622, - 0.15037854, - 1.3167083, - 1.24502001, - 0.42421749, - -0.31524216, - 0.0315209, - 0.52553104, - -0.12301504, - -1.01849236, - 1.66776635, - 0.98537117, - -0.52051714, - 1.15968179, - -0.41569573, - -0.30220252, - -0.24221614, - 1.01834671, - 0.79537263, - 0.47201206, - 0.53230357, - 0.95621545, - 0.67045256, - 0.56275657, - 0.66782958, - 1.04360486, - 1.48513602, - 1.3954355, - 0.72727048, - 0.65556693, - -0.80579428, - -0.6951795, - -0.41933972, - 0.79093036, - -0.13152736, - 0.19672923, - 0.61027041, - 0.34030084, - 1.78314475, - 0.99044934, - 0.51198086, - -0.2394725, - -0.50457759, - 0.57011011, - 1.47181142, - 0.73780165, - 1.67648809, - 0.48628167, - 1.12112297, - 0.49702466, - 0.67702156, - 1.09022274, - 1.06579346, - 0.80085527, - 0.11645159, - -1.58625901, - 1.60302073, - 0.85477186, - -0.06355492, - 0.94015912, - 0.81545671, - 1.36875334, - -0.09772053, - -0.39652235, - 0.11642491, - -0.13081192, - 0.38424794, - 0.07590825, - 0.04093988, - 1.41189211, - 1.21406156, - 0.73694675, - 0.0660595, - 0.25210641, - 0.06578272, - 0.26790674, - 0.77504874, - -0.7038388, - -0.32709336, - 1.25433616, - -0.01871021, - -0.064151, - 0.7459666, - -0.28951973, - -0.8878012, - -0.43049727, - -0.30480856, - 2.2053675, - 0.1190498, - 0.21880221, - 0.77411644, - 1.22875531, - 0.63578569, - 0.96799128, - 0.50535637, - -0.68404005, - 0.66134393, - 0.60923846, - -0.39533348, - 0.70696449, - 0.14433537, - 0.06392199, - 0.53761533, - 0.58593561, - 0.45146962, - 1.16348982, - 0.9917836, - -0.35454457, - 0.52929041, - 0.83350277, + 0.53154128, + 1.69014751, + 0.38863621, + 1.36550742, + 0.54937705, + 0.85452502, + 1.34000193, + 0.04963276, + 0.73207585, + -0.45504452, + 0.99087969, + 0.53800418, + 1.69481083, + 0.19015456, + 0.54587521, + 0.51198155, + 0.17514563, + -0.62023686, + 0.73211032, + 1.0269751, + 0.82004582, + 1.07714066, + 1.27243655, + 0.8603388, + -0.96709536, + -0.4963755, + 0.47436677, + 0.34392296, + 0.08501226, + 0.95779747, + 1.21125461, + -0.04609757, + 0.29714065, + 0.89447118, + 0.00472546, + 0.50365803, + 1.73127064, + 1.04164544, + -0.22236077, + 1.33631993, + 0.96357527, + 1.06122196, + 0.12798557, + 0.4665272, + -0.1162582, + 1.62002463, + 1.44557222, + -0.49218985, + 1.2175313, + 0.25761981, + 0.82879531, + 0.16321047, + 1.34260731, + -0.05709803, + 0.18903618, + 0.76825821, + 0.08211472, + 0.53817434, + 0.53379232, + -0.47094362, + 1.14433914, + 0.03770296, + 1.30737805, + 0.39671022, + 1.22440728, + 0.09600212, + -0.49796137, + -0.44963869, + 0.95208986, + -0.04308567, + 0.45937807, + 2.59887219, + 0.36326674, + 1.16659465, + 2.26888158, + -0.64081701, + 0.13085995, + 1.5847621, + 0.29342994, + -0.7802778, + 0.62631291, + 0.56155063, + 0.63017152, + 1.88801376, + 0.32864795, + 0.19722366, + 0.62506725, + -0.04154236, + 0.74923865, + 0.64958051, + 1.05205509, + 1.12818507, + 0.35463018, + 1.49394093, + -1.32280176, + 0.48922758, + -0.67185953, + 0.01282045, + -0.00832875, + 0.60746178, + 1.04869967, + 0.43197615, + 0.14665959, + -0.08117829, + 0.43216574, + 0.87241428, + -0.07985268, + -0.93380876, + 1.73662159, + 0.23926283, + -0.69068641, + 1.17829179, + -0.16332134, + -0.5112194, + -0.43442261, + 0.34948852, + 1.11002685, + 0.42445302, + 0.68379355, + -0.12877628, + 0.59561974, + 0.67230016, + 1.67895815, + 1.51053172, + 1.14415702, + 1.00682996, + 1.09882483, + 0.28820149, + -0.75250142, + -0.66453929, + -0.0991988, + 0.2907921, + 0.04077525, + -0.44036405, + 0.44894708, + 0.68646345, + 0.03986746, + 0.50061203, + 1.18904715, + 0.36231722, + -0.16347099, + 0.35343108, + 1.15870795, + 0.5973369, + 1.50731862, + 0.69983246, + 1.50854043, + 0.97489667, + 0.25267479, + 0.26369507, + 1.59775053, + 1.56383915, + 0.1721522, + -0.96935772, + 1.47191825, + 0.79858327, + 0.69071774, + -0.17667758, + 0.61438524, + 0.99424152, + -0.23558854, + -0.27873225, + 0.16615446, + 0.02589567, + -0.38007309, + 0.24960815, + 1.17127086, + 1.96577002, + 0.83224965, + 0.9386008, + -0.16018964, + 0.25239747, + -0.09936852, + -0.20376629, + 1.39291948, + -0.2083358, + 0.51435573, + 1.38304537, + 0.23272616, + -0.15257289, + 0.77293718, + 0.33558962, + -0.99534345, + -0.03472376, + 0.07169895, + 1.62726823, + 0.08074445, + 0.38765492, + 0.7844363, + 0.89340893, + 0.28605836, + 0.83632054, + 0.54210362, + -0.55168005, + 0.91756515, + 0.16982656, + -0.36404392, + 1.0011945, + -0.2659181, + 0.31691263, ] ), } From ffed24c613d530d45f9626ec40b1a12d1d535b88 Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Tue, 27 Nov 2018 17:52:42 +0100 Subject: [PATCH 3/4] Removed leftover print statement --- pymc3/distributions/distribution.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index e5681ea34c..0dcb6bdb48 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -512,10 +512,6 @@ def _draw_value(param, point=None, givens=None, size=None): else: output = func(*values) return output - print(param, - type(param), - isinstance(param, tt.TensorVariable), - isinstance(param, (tt.TensorVariable, MultiObservedRV))) raise ValueError('Unexpected type in draw_value: %s' % type(param)) From edd13f7a8acd666f47382e27ab667da8f693cffa Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Thu, 29 Nov 2018 11:05:51 +0100 Subject: [PATCH 4/4] Added release notes and draw values context managers to mixture and multivariate distributions that make many calls to draw_values or other distributions random methods within their own random. --- RELEASE-NOTES.md | 11 +++++++++++ pymc3/distributions/mixture.py | 26 +++++++++++++++----------- pymc3/distributions/multivariate.py | 26 ++++++++++++++------------ 3 files changed, 40 insertions(+), 23 deletions(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index a493e3135a..45b9c3e75d 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -10,12 +10,23 @@ - Add log CDF functions to continuous distributions: `Beta`, `Cauchy`, `ExGaussian`, `Exponential`, `Flat`, `Gumbel`, `HalfCauchy`, `HalfFlat`, `HalfNormal`, `Laplace`, `Logistic`, `Lognormal`, `Normal`, `Pareto`, `StudentT`, `Triangular`, `Uniform`, `Wald`, `Weibull`. - Behavior of `sample_posterior_predictive` is now to produce posterior predictive samples, in order, from all values of the `trace`. Previously, by default it would produce 1 chain worth of samples, using a random selection from the `trace` (#3212) - Show diagnostics for initial energy errors in HMC and NUTS. +- PR #3273 has added the `distributions.distribution._DrawValuesContext` context + manager. This is used to store the values already drawn in nested `random` + and `draw_values` calls, enabling `draw_values` to draw samples from the + joint probability distribution of RVs and not the marginals. Custom + distributions that must call `draw_values` several times in their `random` + method, or that invoke many calls to other distribution's `random` methods + (e.g. mixtures) must do all of these calls under the same `_DrawValuesContext` + context manager instance. If they do not, the conditional relations between + the distribution's parameters could be broken, and `random` could return + values drawn from an incorrect distribution. ### Maintenance - Fixed Triangular distribution `c` attribute handling in `random` and updated sample codes for consistency (#3225) - Renamed `sample_ppc()` and `sample_ppc_w()` to `sample_posterior_predictive()` and `sample_posterior_predictive_w()`, respectively. - Refactor SMC and properly compute marginal likelihood (#3124) +- Fix for #3210. Now `distribution.draw_values(params)`, will draw the `params` values from their joint probability distribution and not from combinations of their marginals (Refer to PR #3273). ## PyMC 3.5 (July 21 2018) diff --git a/pymc3/distributions/mixture.py b/pymc3/distributions/mixture.py index 38f34d6c0a..08a3379be7 100644 --- a/pymc3/distributions/mixture.py +++ b/pymc3/distributions/mixture.py @@ -4,7 +4,8 @@ from pymc3.util import get_variable_name from ..math import logsumexp from .dist_math import bound, random_choice -from .distribution import Discrete, Distribution, draw_values, generate_samples +from .distribution import (Discrete, Distribution, draw_values, + generate_samples, _DrawValuesContext) from .continuous import get_tau_sd, Normal @@ -147,8 +148,9 @@ def logp(self, value): broadcast_conditions=False) def random(self, point=None, size=None): - w = draw_values([self.w], point=point)[0] - comp_tmp = self._comp_samples(point=point, size=None) + with _DrawValuesContext() as draw_context: + w = draw_values([self.w], point=point)[0] + comp_tmp = self._comp_samples(point=point, size=None) if np.asarray(self.shape).size == 0: distshape = np.asarray(np.broadcast(w, comp_tmp).shape)[..., :-1] else: @@ -163,7 +165,8 @@ def random(self, point=None, size=None): dist_shape=distshape, size=size).squeeze() if (size is None) or (distshape.size == 0): - comp_samples = self._comp_samples(point=point, size=size) + with draw_context: + comp_samples = self._comp_samples(point=point, size=size) if comp_samples.ndim > 1: samples = np.squeeze(comp_samples[np.arange(w_samples.size), ..., w_samples]) else: @@ -172,13 +175,14 @@ def random(self, point=None, size=None): if w_samples.ndim == 1: w_samples = np.reshape(np.tile(w_samples, size), (size,) + w_samples.shape) samples = np.zeros((size,)+tuple(distshape)) - for i in range(size): - w_tmp = w_samples[i, :] - comp_tmp = self._comp_samples(point=point, size=None) - if comp_tmp.ndim > 1: - samples[i, :] = np.squeeze(comp_tmp[np.arange(w_tmp.size), ..., w_tmp]) - else: - samples[i, :] = np.squeeze(comp_tmp[w_tmp]) + with draw_context: + for i in range(size): + w_tmp = w_samples[i, :] + comp_tmp = self._comp_samples(point=point, size=None) + if comp_tmp.ndim > 1: + samples[i, :] = np.squeeze(comp_tmp[np.arange(w_tmp.size), ..., w_tmp]) + else: + samples[i, :] = np.squeeze(comp_tmp[w_tmp]) return samples diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 1880bb2996..f50a25e287 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -16,7 +16,8 @@ from pymc3.theanof import floatX from . import transforms from pymc3.util import get_variable_name -from .distribution import Continuous, Discrete, draw_values, generate_samples +from .distribution import (Continuous, Discrete, draw_values, generate_samples, + _DrawValuesContext) from ..model import Deterministic from .continuous import ChiSquared, Normal from .special import gammaln, multigammaln @@ -338,18 +339,19 @@ def __init__(self, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, self.mean = self.median = self.mode = self.mu = self.mu def random(self, point=None, size=None): - nu, mu = draw_values([self.nu, self.mu], point=point, size=size) - if self._cov_type == 'cov': - cov, = draw_values([self.cov], point=point, size=size) - dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov) - elif self._cov_type == 'tau': - tau, = draw_values([self.tau], point=point, size=size) - dist = MvNormal.dist(mu=np.zeros_like(mu), tau=tau) - else: - chol, = draw_values([self.chol_cov], point=point, size=size) - dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol) + with _DrawValuesContext(): + nu, mu = draw_values([self.nu, self.mu], point=point, size=size) + if self._cov_type == 'cov': + cov, = draw_values([self.cov], point=point, size=size) + dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov) + elif self._cov_type == 'tau': + tau, = draw_values([self.tau], point=point, size=size) + dist = MvNormal.dist(mu=np.zeros_like(mu), tau=tau) + else: + chol, = draw_values([self.chol_cov], point=point, size=size) + dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol) - samples = dist.random(point, size) + samples = dist.random(point, size) chi2 = np.random.chisquare return (np.sqrt(nu) * samples.T / chi2(nu, size)).T + mu