From e99207ddf95ec0e7b6b9733711086738f36bb3de Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Mon, 26 Feb 2018 17:54:29 -0500 Subject: [PATCH 1/9] Add sample_prior function --- pymc3/sampling.py | 59 +++++++++++++++++++++++++++++++++++- pymc3/tests/test_sampling.py | 12 ++++++++ 2 files changed, 70 insertions(+), 1 deletion(-) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 067c17a890..1ca0a02d86 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -25,7 +25,7 @@ import sys sys.setrecursionlimit(10000) -__all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts'] +__all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts', 'sample_prior'] STEP_METHODS = (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, Slice, CategoricalGibbsMetropolis) @@ -1206,6 +1206,63 @@ def sample_ppc_w(traces, samples=None, models=None, weights=None, return {k: np.asarray(v) for k, v in ppc.items()} +def sample_prior(samples=500, model=None, vars=None, size=None, + random_seed=None, progressbar=True): + """Generate samples from the prior of a model. + + Parameters + ---------- + samples : int + Number of samples from the prior to generate. Defaults to 500. + model : Model (optional if in `with` context) + vars : iterable + Variables for which to compute the posterior predictive samples. + Defaults to `model.named_vars`. + size : int + The number of random draws from the distribution specified by the + parameters in each sample of the trace. + random_seed : int + Seed for the random number generator. + progressbar : bool + Whether or not to display a progress bar in the command line. + + Returns + ------- + dict + Dictionary with the variables as keys. The values are arrays of prior samples. + """ + + model = modelcontext(model) + + if vars is None: + vars = set(model.named_vars.keys()) + + if random_seed is not None: + np.random.seed(random_seed) + + if progressbar: + indices = tqdm(range(samples)) + + try: + prior = {var: [] for var in vars} + for _ in indices: + point = {} + for var_name, var in model.named_vars.items(): + val = var.distribution.random(point=point, size=size) + point[var_name] = val + if var_name in vars: + prior[var_name].append(val) + + except KeyboardInterrupt: + pass + + finally: + if progressbar: + indices.close() + + return {k: np.asarray(v) for k, v in prior.items()} + + def init_nuts(init='auto', chains=1, n_init=500000, model=None, random_seed=None, progressbar=True, **kwargs): """Set up the mass matrix initialization for NUTS. diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py index 30362019c0..f3810c72b8 100644 --- a/pymc3/tests/test_sampling.py +++ b/pymc3/tests/test_sampling.py @@ -302,3 +302,15 @@ def test_exec_nuts_init(method): assert len(start) == 2 assert isinstance(start[0], dict) assert 'a' in start[0] and 'b_log__' in start[0] + + +def test_sample_prior(): + observed = np.random.normal(10, 1, size=200) + with pm.Model(): + # Use a prior that's way off to show we're actually sampling from it + mu = pm.Normal('mu', mu=-10, sd=1) + pm.Normal('x_obs', mu=mu, sd=1, observed=observed) + prior = pm.sample_prior() + + assert (prior['mu'] < 0).all() + assert (prior['x_obs'] < 0).all() \ No newline at end of file From 6e75ed6774d21f9ec8cd542f0e30b30ff8902f58 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Tue, 27 Feb 2018 21:38:42 -0500 Subject: [PATCH 2/9] Allow sample_prior with transformed distributions --- pymc3/sampling.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 1ca0a02d86..f8f55876e1 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -15,7 +15,7 @@ from .step_methods import (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, Slice, CompoundStep, arraystep) -from .util import update_start_vals +from .util import update_start_vals, is_transformed_name from .vartypes import discrete_types from pymc3.step_methods.hmc import quadpotential from pymc3 import plots @@ -1248,7 +1248,11 @@ def sample_prior(samples=500, model=None, vars=None, size=None, for _ in indices: point = {} for var_name, var in model.named_vars.items(): - val = var.distribution.random(point=point, size=size) + if is_transformed_name(var_name): + trans = var.distribution.transform_used.forward + val = trans(var.distribution.dist.random(point=point, size=size)).eval() + else: + val = var.distribution.random(point=point, size=size) point[var_name] = val if var_name in vars: prior[var_name].append(val) From 99f71a3c67adc678dba837402d92c4fd5ec18cbb Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Tue, 27 Feb 2018 21:38:55 -0500 Subject: [PATCH 3/9] Update and tidy release notes --- RELEASE-NOTES.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index cb27ab8538..78792e375e 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -11,8 +11,7 @@ Gaussian Process implementation for efficient inference, along with lower-level functions such as `cartesian` and `kronecker` products. - Added `Coregion` covariance function. -- Add new 'pairplot' function, for plotting scatter or hexbin matrices of sampled parameters. - Optionally it can plot divergences. +- Add new 'pairplot' function, for plotting scatter or hexbin matrices of sampled parameters. Optionally it can plot divergences. - Plots of discrete distributions in the docstrings - Add logitnormal distribution - Densityplot: add support for discrete variables @@ -21,6 +20,7 @@ - Changed the `compare` function to accept a dictionary of model-trace pairs instead of two separate lists of models and traces. - add test and support for creating multivariate mixture and mixture of mixtures - `distribution.draw_values`, now is also able to draw values from conditionally dependent RVs, such as autotransformed RVs (Refer to PR #2902). +- New function `pm.sample_prior` which generates test data from a model in the absence of data. ### Fixes From 8c0d9b3eac8178e03c3238e9b7b6e8fe556ab36f Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Wed, 28 Feb 2018 08:59:03 -0500 Subject: [PATCH 4/9] Fix speed problem --- pymc3/sampling.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index f8f55876e1..0430c493c8 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -1249,8 +1249,7 @@ def sample_prior(samples=500, model=None, vars=None, size=None, point = {} for var_name, var in model.named_vars.items(): if is_transformed_name(var_name): - trans = var.distribution.transform_used.forward - val = trans(var.distribution.dist.random(point=point, size=size)).eval() + val = var.distribution.dist.random(point=point, size=size) else: val = var.distribution.random(point=point, size=size) point[var_name] = val From ea3026bd8f45f6bf7519a8418934804fc93e0764 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Sun, 4 Mar 2018 09:42:27 -0500 Subject: [PATCH 5/9] Works with deterministic variables now --- pymc3/sampling.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 0430c493c8..56aac38b31 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -1248,10 +1248,14 @@ def sample_prior(samples=500, model=None, vars=None, size=None, for _ in indices: point = {} for var_name, var in model.named_vars.items(): - if is_transformed_name(var_name): - val = var.distribution.dist.random(point=point, size=size) + if hasattr(var, 'distribution'): + if is_transformed_name(var_name): + val = var.distribution.dist.random(point=point, size=size) + else: + val = var.distribution.random(point=point, size=size) else: - val = var.distribution.random(point=point, size=size) + val = var.eval({model.named_vars[v]: point[v] for v in pm.model.get_named_nodes(var)}) + point[var_name] = val if var_name in vars: prior[var_name].append(val) From 16c0b3ab8a8d35ab4b6ff1fbeb09883e5b23b084 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Sun, 4 Mar 2018 09:45:04 -0500 Subject: [PATCH 6/9] Update test --- pymc3/tests/test_sampling.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py index f3810c72b8..b897ee782b 100644 --- a/pymc3/tests/test_sampling.py +++ b/pymc3/tests/test_sampling.py @@ -309,8 +309,10 @@ def test_sample_prior(): with pm.Model(): # Use a prior that's way off to show we're actually sampling from it mu = pm.Normal('mu', mu=-10, sd=1) + positive_mu = pm.Deterministic('positive_mu', np.abs(mu)) pm.Normal('x_obs', mu=mu, sd=1, observed=observed) prior = pm.sample_prior() assert (prior['mu'] < 0).all() + assert (prior['positive_mu'] > 0).all() assert (prior['x_obs'] < 0).all() \ No newline at end of file From 9a80d0d4cd5c567602432cdaa393447d8fc69a31 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Sun, 25 Mar 2018 15:45:39 -0400 Subject: [PATCH 7/9] Update sampling to use `draw_values` --- pymc3/distributions/distribution.py | 45 ++++++++++++++++------ pymc3/sampling.py | 59 ++++++++++------------------- pymc3/tests/test_sampling.py | 4 +- 3 files changed, 56 insertions(+), 52 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 8139c8dd46..a800b534a5 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -210,7 +210,7 @@ def random(self, *args, **kwargs): -def draw_values(params, point=None): +def draw_values(params, point=None, size=None): """ Draw (fix) parameter values. Handles a number of cases: @@ -254,7 +254,7 @@ def draw_values(params, point=None): # Init givens and the stack of nodes to try to `_draw_value` from givens = {} - stored = set([]) # Some nodes + stored = set() # Some nodes stack = list(leaf_nodes.values()) # A queue would be more appropriate while stack: next_ = stack.pop(0) @@ -279,13 +279,14 @@ def draw_values(params, point=None): # 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.keys() if k in children] + 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)) + 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 @@ -295,10 +296,27 @@ def draw_values(params, point=None): if node is not None and node.name not in stored and node not in params]) - values = [] - for param in params: - values.append(_draw_value(param, point=point, givens=givens.values())) - return values + + # Funny dance to resolve missing nodes + 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(j) for j in to_eval])) + to_eval = set(missing_inputs) + missing_inputs = set() + for param in to_eval: + try: # might evaluate in a bad order, + evaluated[param] = _draw_value(params[param], point=point, givens=givens.values(), size=size) + if any(param in j for j in named_nodes_children.values()): + givens[param.name] = (params[param], evaluated[param]) + except theano.gof.fg.MissingInputError: + missing_inputs.add(param) + + + return [evaluated[j] for j in params] # set the order back @memoize @@ -326,7 +344,7 @@ def _compile_theano_function(param, vars, givens=None): allow_input_downcast=True) -def _draw_value(param, point=None, givens=None): +def _draw_value(param, point=None, givens=None, size=None): """Draw a random value from a distribution or return a constant. Parameters @@ -355,14 +373,19 @@ def _draw_value(param, point=None, givens=None): if point and hasattr(param, 'model') and param.name in point: return point[param.name] elif hasattr(param, 'random') and param.random is not None: - return param.random(point=point, size=None) + return param.random(point=point, size=size) + elif hasattr(param, 'distribution') and hasattr(param.distribution, 'random') and param.distribution.random is not None: + return param.distribution.random(point=point, size=size) else: if givens: variables, values = list(zip(*givens)) else: variables = values = [] func = _compile_theano_function(param, variables) - return func(*values) + if size is None: + return func(*values) + else: + return np.array([func(*value) for value in zip(*values)]) else: raise ValueError('Unexpected type in draw_value: %s' % type(param)) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 56aac38b31..38bf66593c 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -11,11 +11,12 @@ from .backends.base import BaseTrace, MultiTrace from .backends.ndarray import NDArray +from .distributions import draw_values from .model import modelcontext, Point, all_continuous from .step_methods import (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, Slice, CompoundStep, arraystep) -from .util import update_start_vals, is_transformed_name +from .util import update_start_vals, is_transformed_name, get_untransformed_name, get_default_varnames from .vartypes import discrete_types from pymc3.step_methods.hmc import quadpotential from pymc3 import plots @@ -25,7 +26,7 @@ import sys sys.setrecursionlimit(10000) -__all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts', 'sample_prior'] +__all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts', 'sample_generative'] STEP_METHODS = (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, Slice, CategoricalGibbsMetropolis) @@ -1206,9 +1207,8 @@ def sample_ppc_w(traces, samples=None, models=None, weights=None, return {k: np.asarray(v) for k, v in ppc.items()} -def sample_prior(samples=500, model=None, vars=None, size=None, - random_seed=None, progressbar=True): - """Generate samples from the prior of a model. +def sample_generative(samples=500, model=None, vars=None, random_seed=None): + """Generate samples from the generative model. Parameters ---------- @@ -1218,20 +1218,14 @@ def sample_prior(samples=500, model=None, vars=None, size=None, vars : iterable Variables for which to compute the posterior predictive samples. Defaults to `model.named_vars`. - size : int - The number of random draws from the distribution specified by the - parameters in each sample of the trace. random_seed : int Seed for the random number generator. - progressbar : bool - Whether or not to display a progress bar in the command line. Returns ------- dict Dictionary with the variables as keys. The values are arrays of prior samples. """ - model = modelcontext(model) if vars is None: @@ -1240,34 +1234,21 @@ def sample_prior(samples=500, model=None, vars=None, size=None, if random_seed is not None: np.random.seed(random_seed) - if progressbar: - indices = tqdm(range(samples)) - - try: - prior = {var: [] for var in vars} - for _ in indices: - point = {} - for var_name, var in model.named_vars.items(): - if hasattr(var, 'distribution'): - if is_transformed_name(var_name): - val = var.distribution.dist.random(point=point, size=size) - else: - val = var.distribution.random(point=point, size=size) - else: - val = var.eval({model.named_vars[v]: point[v] for v in pm.model.get_named_nodes(var)}) - - point[var_name] = val - if var_name in vars: - prior[var_name].append(val) - - except KeyboardInterrupt: - pass - - finally: - if progressbar: - indices.close() - - return {k: np.asarray(v) for k, v in prior.items()} + names = get_default_varnames(model.named_vars, include_transformed=False) + # draw_values fails with auto-transformed variables. transform them later! + values = draw_values([model[name] for name in names], size=samples) + + data = {k: v for k, v in zip(names, values)} + + prior = {} + for var_name in vars: + if var_name in data: + prior[var_name] = data[var_name] + elif is_transformed_name(var_name): + untransformed = get_untransformed_name(var_name) + if untransformed in data: + prior[var_name] = model[untransformed].transformation.forward(data[untransformed]).eval() + return prior def init_nuts(init='auto', chains=1, n_init=500000, model=None, diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py index b897ee782b..6b4d639246 100644 --- a/pymc3/tests/test_sampling.py +++ b/pymc3/tests/test_sampling.py @@ -304,14 +304,14 @@ def test_exec_nuts_init(method): assert 'a' in start[0] and 'b_log__' in start[0] -def test_sample_prior(): +def test_sample_generative(): observed = np.random.normal(10, 1, size=200) with pm.Model(): # Use a prior that's way off to show we're actually sampling from it mu = pm.Normal('mu', mu=-10, sd=1) positive_mu = pm.Deterministic('positive_mu', np.abs(mu)) pm.Normal('x_obs', mu=mu, sd=1, observed=observed) - prior = pm.sample_prior() + prior = pm.sample_generative() assert (prior['mu'] < 0).all() assert (prior['positive_mu'] > 0).all() From b45d6e7d0b577bf840676d7b37fb588456f8c680 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Sun, 25 Mar 2018 23:13:55 -0400 Subject: [PATCH 8/9] Actually add givens --- pymc3/distributions/distribution.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index a800b534a5..67e079d170 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -304,14 +304,14 @@ def draw_values(params, point=None, size=None): missing_inputs = set(params) while to_eval or missing_inputs: if to_eval == missing_inputs: - raise ValueError('Cannot resolve inputs for {}'.format([str(j) for j in to_eval])) + 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 in to_eval: try: # might evaluate in a bad order, evaluated[param] = _draw_value(params[param], point=point, givens=givens.values(), size=size) - if any(param in j for j in named_nodes_children.values()): - givens[param.name] = (params[param], evaluated[param]) + if any(params[param] in j for j in named_nodes_children.values()): + givens[params[param].name] = (params[param], evaluated[param]) except theano.gof.fg.MissingInputError: missing_inputs.add(param) From 4541d16dad0ff9a85d7664f3742b6262cee51ee4 Mon Sep 17 00:00:00 2001 From: Colin Carroll Date: Mon, 26 Mar 2018 11:37:26 -0400 Subject: [PATCH 9/9] Incorporate comments, make more readable --- pymc3/distributions/distribution.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 67e079d170..68e70cf017 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -1,3 +1,4 @@ +import collections import numbers import numpy as np import theano.tensor as tt @@ -307,13 +308,14 @@ def draw_values(params, point=None, size=None): 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 in to_eval: + for param_idx in to_eval: + param = params[param_idx] try: # might evaluate in a bad order, - evaluated[param] = _draw_value(params[param], point=point, givens=givens.values(), size=size) - if any(params[param] in j for j in named_nodes_children.values()): - givens[params[param].name] = (params[param], evaluated[param]) + 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]) except theano.gof.fg.MissingInputError: - missing_inputs.add(param) + missing_inputs.add(param_idx) return [evaluated[j] for j in params] # set the order back