diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 69c44cd3b82..3572ada903e 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -5,6 +5,9 @@ - Improve NUTS initialization `advi+adapt_diag_grad` and add `jitter+adapt_diag_grad` (#2643) - Fixed `compareplot` to use `loo` output. +### New Features +- Michael Osthege added support for population-samplers and implemented differential evolution metropolis (`DEMetropolis`). For models with correlated dimensions that can not use gradient-based samplers, the `DEMetropolis` sampler can give higher effective sampling rates. (also see [PR#2735](https://github.com/pymc-devs/pymc3/pull/2735)) + ## PyMC3 3.2 (October 10, 2017) diff --git a/pymc3/examples/samplers_mvnormal.py b/pymc3/examples/samplers_mvnormal.py new file mode 100644 index 00000000000..14b87d7958c --- /dev/null +++ b/pymc3/examples/samplers_mvnormal.py @@ -0,0 +1,100 @@ +""" +Comparing different samplers on a correlated bivariate normal distribution. + +This example will sample a bivariate normal with Metropolis, NUTS and DEMetropolis +at two correlations (0, 0.9) and print out the effective sample sizes, runtime and +normalized effective sampling rates. +""" + + +import numpy as np +import time +import pandas as pd +import pymc3 as pm +import theano.tensor as tt + +# with this flag one can switch between defining the bivariate normal as +# either a 2D MvNormal (USE_XY = False) split up the two dimensions into +# two variables 'x' and 'y'. The latter is recommended because it highlights +# different behaviour with respect to blocking. +USE_XY = True + +def run(steppers, p): + steppers = set(steppers) + traces = {} + effn = {} + runtimes = {} + + with pm.Model() as model: + if USE_XY: + x = pm.Flat('x') + y = pm.Flat('y') + mu = np.array([0.,0.]) + cov = np.array([[1.,p],[p,1.]]) + z = pm.MvNormal.dist(mu=mu, cov=cov, shape=(2,)).logp(tt.stack([x,y])) + pot = pm.Potential('logp_xy', z) + start = {'x': 0, 'y': 0} + else: + mu = np.array([0.,0.]) + cov = np.array([[1.,p],[p,1.]]) + z = pm.MvNormal('z', mu=mu, cov=cov, shape=(2,)) + start={'z': [0, 0]} + + for step_cls in steppers: + name = step_cls.__name__ + t_start = time.time() + mt = pm.sample( + draws=10000, + chains=16, parallelize=False, + step=step_cls(), + start=start + ) + runtimes[name] = time.time() - t_start + print('{} samples across {} chains'.format(len(mt) * mt.nchains, mt.nchains)) + traces[name] = mt + en = pm.diagnostics.effective_n(mt) + print('effective: {}\r\n'.format(en)) + if USE_XY: + effn[name] = np.mean(en['x']) / len(mt) / mt.nchains + else: + effn[name] = np.mean(en['z']) / len(mt) / mt.nchains + return traces, effn, runtimes + + +if __name__ == '__main__': + methods = [ + pm.Metropolis, + pm.Slice, + pm.NUTS, + pm.DEMetropolis + ] + names = [c.__name__ for c in methods] + + df_base = pd.DataFrame(columns=['p'] + names) + df_base['p'] = [.0,.9] + df_base = df_base.set_index('p') + + df_effectiven = df_base.copy() + df_runtime = df_base.copy() + df_performance = df_base.copy() + + for p in df_effectiven.index: + trace, rate, runtime = run(methods, p) + for name in names: + df_effectiven.set_value(p, name, rate[name]) + df_runtime.set_value(p, name, runtime[name]) + df_performance.set_value(p, name, rate[name] / runtime[name]) + + print('\r\nEffective sample size [0...1]') + print(df_effectiven.T.to_string(float_format='{:.3f}'.format)) + + print('\r\nRuntime [s]') + print(df_runtime.T.to_string(float_format='{:.1f}'.format)) + + if 'NUTS' in names: + print('\r\nNormalized effective sampling rate [0...1]') + df_performance = df_performance.T / df_performance.loc[0]['NUTS'] + else: + print('\r\nNormalized effective sampling rate [1/s]') + df_performance = df_performance.T + print(df_performance.to_string(float_format='{:.3f}'.format)) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 16db40b9515..f2bd3865fa8 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -1,5 +1,7 @@ from collections import defaultdict, Iterable +from copy import copy import pickle +from six import integer_types from joblib import Parallel, delayed import numpy as np @@ -10,9 +12,9 @@ from .backends.base import BaseTrace, MultiTrace from .backends.ndarray import NDArray from .model import modelcontext, Point -from .step_methods import (NUTS, HamiltonianMC, SGFS, Metropolis, BinaryMetropolis, +from .step_methods import (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, - Slice, CompoundStep) + Slice, CompoundStep, arraystep) from .util import update_start_vals from .vartypes import discrete_types from pymc3.step_methods.hmc import quadpotential @@ -23,7 +25,7 @@ __all__ = ['sample', 'iter_sample', 'sample_ppc', 'sample_ppc_w', 'init_nuts'] -STEP_METHODS = (NUTS, HamiltonianMC, SGFS, Metropolis, BinaryMetropolis, +STEP_METHODS = (NUTS, HamiltonianMC, Metropolis, BinaryMetropolis, BinaryGibbsMetropolis, Slice, CategoricalGibbsMetropolis) @@ -143,6 +145,19 @@ def assign_step_methods(model, step=None, methods=STEP_METHODS, return instantiate_steppers(model, steps, selected_steps, step_kwargs) +def print_step_hierarchy(s, level=0): + if isinstance(s, (list, tuple)): + pm._log.info('>' * level + 'list') + for i in s: + print_step_hierarchy(i, level+1) + elif isinstance(s, CompoundStep): + pm._log.info('>' * level + 'CompoundStep') + for i in s.methods: + print_step_hierarchy(i, level+1) + else: + pm._log.info('>' * level + '{}: {}'.format(s.__class__.__name__, s.vars)) + + def _cpu_count(): """Try to guess the number of CPUs in the system. @@ -357,8 +372,10 @@ def sample(draws=500, step=None, init='auto', n_init=200000, start=None, else: step = assign_step_methods(model, step, step_kwargs=step_kwargs) + if isinstance(step, list): + step = CompoundStep(step) if start is None: - start = [None] * chains + start = {} if isinstance(start, dict): start = [start] * chains @@ -380,23 +397,36 @@ def sample(draws=500, step=None, init='auto', n_init=200000, start=None, sample_args.update(kwargs) - parallel = njobs > 1 and chains > 1 + has_population_samplers = np.any([ + isinstance(m, arraystep.PopulationArrayStepShared) + for m in (step.methods if isinstance(step, CompoundStep) else [step]) + ]) + parallel = njobs > 1 and chains > 1 and not has_population_samplers if parallel: + pm._log.info('Multiprocess sampling ({} chains in {} jobs)'.format(chains, njobs)) + print_step_hierarchy(step) try: trace = _mp_sample(**sample_args) except pickle.PickleError: - pm._log.warn("Could not pickle model, sampling sequentially.") + pm._log.warn("Could not pickle model, sampling singlethreaded.") pm._log.debug('Pickling error:', exec_info=True) parallel = False except AttributeError as e: if str(e).startswith("AttributeError: Can't pickle"): - pm._log.warn("Could not pickle model, sampling sequentially.") + pm._log.warn("Could not pickle model, sampling singlethreaded.") pm._log.debug('Pickling error:', exec_info=True) parallel = False else: raise if not parallel: - trace = _sample_many(**sample_args) + if has_population_samplers: + pm._log.info('Population sampling ({} chains)'.format(chains)) + print_step_hierarchy(step) + trace = _sample_population(**sample_args) + else: + pm._log.info('Sequential sampling ({} chains in 1 job)'.format(chains)) + print_step_hierarchy(step) + trace = _sample_many(**sample_args) discard = tune if discard_tuned_samples else 0 return trace[discard:] @@ -448,6 +478,23 @@ def _sample_many(draws, chain, chains, start, random_seed, **kwargs): return MultiTrace(traces) +def _sample_population(draws, chain, chains, start, random_seed, step, tune, + model, progressbar=None, parallelize=False, **kwargs): + # create the generator that iterates all chains in parallel + chains = [chain + c for c in range(chains)] + sampling = _prepare_iter_population(draws, chains, step, start, parallelize, + tune=tune, model=model, random_seed=random_seed) + + if progressbar: + sampling = tqdm(sampling, total=draws) + + latest_traces = None + for it,traces in enumerate(sampling): + latest_traces = traces + # TODO: add support for liveplot during population-sampling + return MultiTrace(latest_traces) + + def _sample(chain, progressbar, random_seed, start, draws=None, step=None, trace=None, tune=None, model=None, live_plot=False, live_plot_kwargs=None, **kwargs): @@ -580,6 +627,267 @@ def _iter_sample(draws, step, start=None, trace=None, chain=0, tune=None, step.report._finalize(strace) +class PopulationStepper(object): + def __init__(self, steppers, parallelize): + """Tries to use multiprocessing to parallelize chains. + + Falls back to sequential evaluation if multiprocessing fails. + + In the multiprocessing mode of operation, a new process is started for each + chain/stepper and Pipes are used to communicate with the main process. + + Parameters + ---------- + steppers : list + A collection of independent step methods, one for each chain. + parallelize : bool + Indicates if chain parallelization is desired + """ + self.nchains = len(steppers) + self.is_parallelized = False + self._master_ends = [] + self._processes = [] + self._steppers = steppers + if parallelize and sys.version_info >= (3,4): + try: + # configure a child process for each stepper + pm._log.info('Attempting to parallelize chains.') + import multiprocessing + for c,stepper in enumerate(steppers): + slave_end, master_end = multiprocessing.Pipe() + stepper_dumps = pickle.dumps(stepper, protocol=4) + process = multiprocessing.Process( + target=self.__class__._run_slave, + args=(c, stepper_dumps, slave_end), + name='ChainWalker{}'.format(c) + ) + # Starting the process might fail and takes time. + # By doing it in the constructor, the sampling progress bar + # will not be confused by the process start. + process.start() + self._master_ends.append(master_end) + self._processes.append(process) + self.is_parallelized = True + except: + pm._log.info('Population parallelization failed. ' \ + 'Falling back to sequential stepping of chains.') + else: + if parallelize: + warnings.warn('Population parallelization is only supported on Python 3.4 and ' \ + 'higher. All {} chains will step on one process.'.format(self.nchains)) + else: + pm._log.info('Chains are not parallelized. You can enable this by passing ' \ + 'pm.sample(parallelize=True).') + return super(PopulationStepper, self).__init__() + + def __enter__(self): + """Does nothing because processes are already started in __init__.""" + return + + def __exit__(self, exc_type, exc_val, exc_tb): + if len(self._processes) > 0: + try: + for master_end in self._master_ends: + master_end.send(None) + for process in self._processes: + process.join(timeout=3) + except: + pm._log.warning('Termination failed.') + return + + @staticmethod + def _run_slave(c, stepper_dumps, slave_end): + """Started on a separate process to perform stepping of a chain. + + Parameters + ---------- + c : int + number of this chain + stepper : BlockedStep + a step method such as CompoundStep + slave_end : multiprocessing.connection.PipeConnection + This is our connection to the main process + """ + try: + stepper = pickle.loads(stepper_dumps) + # the stepper is not necessarily a PopulationArraySharedStep itself, + # but rather a CompoundStep. PopulationArrayStepShared.population + # has to be updated, therefore we identify the substeppers first. + population_steppers = [] + for sm in (stepper.methods if isinstance(stepper, CompoundStep) else [stepper]): + if isinstance(sm, arraystep.PopulationArrayStepShared): + population_steppers.append(sm) + while True: + incoming = slave_end.recv() + # receiving a None is the signal to exit + if incoming is None: + break + tune_stop, population = incoming + if tune_stop: + stop_tuning(stepper) + # forward the population to the PopulationArrayStepShared objects + # This is necessary because due to the process fork, the population + # object is no longer shared between the steppers. + for popstep in population_steppers: + popstep.population = population + update = stepper.step(population[c]) + slave_end.send(update) + except Exception: + pm._log.exception('ChainWalker{}'.format(c)) + return + + def step(self, tune_stop, population): + """Steps the entire population of chains. + + Parameters + ---------- + tune_stop : bool + Indicates if the condition (i == tune) is fulfilled + population : list + Current Points of all chains + + Returns + ------- + update : Point + The new positions of the chains + """ + updates = [None] * self.nchains + if self.is_parallelized: + for c in range(self.nchains): + self._master_ends[c].send((tune_stop, population)) + # Blockingly get the step outcomes + for c in range(self.nchains): + updates[c] = self._master_ends[c].recv() + else: + for c in range(self.nchains): + if tune_stop: + self._steppers[c] = stop_tuning(self._steppers[c]) + updates[c] = self._steppers[c].step(population[c]) + return updates + + +def _prepare_iter_population(draws, chains, step, start, parallelize, tune=None, + model=None, random_seed=None): + """Prepares a PopulationStepper and traces for population sampling. + + Returns + ------- + _iter_population : generator + The generator the yields traces of all chains at the same time + """ + # chains contains the chain numbers, but for indexing we need indices... + nchains = len(chains) + model = modelcontext(model) + draws = int(draws) + if random_seed is not None: + np.random.seed(random_seed) + if draws < 1: + raise ValueError('Argument `draws` should be above 0.') + + # The initialization of traces, samplers and points must happen in the right order: + # 1. traces are initialized and update_start_vals configures variable transforms + # 2. population of points is created + # 3. steppers are initialized and linked to the points object + # 4. traces are configured to track the sampler stats + # 5. a PopulationStepper is configured for parallelized stepping + + # 1. prepare a BaseTrace for each chain + traces = [_choose_backend(None, chain, model=model) for chain in chains] + for c,strace in enumerate(traces): + # initialize the trace size and variable transforms + if len(strace) > 0: + update_start_vals(start[c], strace.point(-1), model) + else: + update_start_vals(start[c], model.test_point, model) + + # 2. create a population (points) that tracks each chain + # it is updated as the chains are advanced + population = [Point(start[c], model=model) for c in range(nchains)] + + # 3. Set up the steppers + steppers = [None] * nchains + for c in range(nchains): + # need indepenent samplers for each chain + # it is important to copy the actual steppers (but not the delta_logp) + if isinstance(step, CompoundStep): + chainstep = CompoundStep([copy(m) for m in step.methods]) + else: + chainstep = copy(step) + # link population samplers to the shared population state + for sm in (chainstep.methods if isinstance(step, CompoundStep) else [chainstep]): + if isinstance(sm, arraystep.PopulationArrayStepShared): + sm.link_population(population, c) + steppers[c] = chainstep + + # 4. configure tracking of sampler stats + for c in range(nchains): + if steppers[c].generates_stats and traces[c].supports_sampler_stats: + traces[c].setup(draws, c, steppers[c].stats_dtypes) + else: + traces[c].setup(draws, c) + + # 5. configure the PopulationStepper (expensive call) + popstep = PopulationStepper(steppers, parallelize) + + # Because the preperations above are expensive, the actual iterator is + # in another method. This way the progbar will not be disturbed. + return _iter_population(draws, tune, popstep, steppers, traces, population) + + +def _iter_population(draws, tune, popstep, steppers, traces, points): + """Generator that iterates a PopulationStepper. + + Parameters + ---------- + draws : int + number of draws per chain + tune : int + number of tuning steps + popstep : PopulationStepper + the helper object for (parallelized) stepping of chains + steppers : list + The step methods for each chain + traces : list + Traces for each chain + points : list + population of chain states + """ + try: + with popstep: + # iterate draws of all chains + for i in range(draws): + updates = popstep.step(i == tune, points) + + # apply the update to the points and record to the traces + for c,strace in enumerate(traces): + if steppers[c].generates_stats: + points[c], states = updates[c] + if strace.supports_sampler_stats: + strace.record(points[c], states) + else: + strace.record(points[c]) + else: + points[c] = updates[c] + strace.record(points[c]) + # yield the state of all chains in parallel + yield traces + except KeyboardInterrupt: + for c,strace in enumerate(traces): + strace.close() + if hasattr(steppers[c], 'report'): + steppers[c].report._finalize(strace) + raise + except BaseException: + for c,strace in enumerate(traces): + strace.close() + raise + else: + for c,strace in enumerate(traces): + strace.close() + if hasattr(steppers[c], 'report'): + steppers[c].report._finalize(strace) + + def _choose_backend(trace, chain, shortcuts=None, **kwds): if isinstance(trace, BaseTrace): return trace diff --git a/pymc3/step_methods/__init__.py b/pymc3/step_methods/__init__.py index 67212508d3e..d18e36a153c 100644 --- a/pymc3/step_methods/__init__.py +++ b/pymc3/step_methods/__init__.py @@ -3,6 +3,7 @@ from .hmc import HamiltonianMC, NUTS from .metropolis import Metropolis +from .metropolis import DEMetropolis from .metropolis import BinaryMetropolis from .metropolis import BinaryGibbsMetropolis from .metropolis import CategoricalGibbsMetropolis diff --git a/pymc3/step_methods/arraystep.py b/pymc3/step_methods/arraystep.py index 564d37e2fc2..5d6d1db1a25 100644 --- a/pymc3/step_methods/arraystep.py +++ b/pymc3/step_methods/arraystep.py @@ -150,19 +150,58 @@ def __init__(self, vars, shared, blocked=True): self.ordering = ArrayOrdering(vars) self.shared = {str(var): shared for var, shared in shared.items()} self.blocked = blocked + self.bij = None def step(self, point): for var, share in self.shared.items(): share.set_value(point[var]) - bij = DictToArrayBijection(self.ordering, point) + self.bij = DictToArrayBijection(self.ordering, point) if self.generates_stats: - apoint, stats = self.astep(bij.map(point)) - return bij.rmap(apoint), stats + apoint, stats = self.astep(self.bij.map(point)) + return self.bij.rmap(apoint), stats else: - apoint = self.astep(bij.map(point)) - return bij.rmap(apoint) + apoint = self.astep(self.bij.map(point)) + return self.bij.rmap(apoint) + + +class PopulationArrayStepShared(ArrayStepShared): + """Version of ArrayStepShared that allows samplers to access the states + of other chains in the population. + + Works by linking a list of Points that is updated as the chains are iterated. + """ + + def __init__(self, vars, shared, blocked=True): + """ + Parameters + ---------- + vars : list of sampling variables + shared : dict of theano variable -> shared variable + blocked : Boolean (default True) + """ + self.population = None + self.this_chain = None + self.other_chains = None + return super(PopulationArrayStepShared, self).__init__(vars, shared, blocked) + + def link_population(self, population, chain_index): + """Links the sampler to the population. + + Parameters + ---------- + population : list of Points. (The elements of this list must be + replaced with current chain states in every iteration.) + chain_index : int of the index of this sampler in the population + """ + self.population = population + self.this_chain = chain_index + self.other_chains = [c for c in range(len(population)) if c != chain_index] + if not len(self.other_chains) > 1: + raise ValueError('Population is just {} + {}. This is too small. You should ' \ + 'increase the number of chains.'.format(self.this_chain, self.other_chains)) + return class GradientSharedStep(BlockedStep): diff --git a/pymc3/step_methods/metropolis.py b/pymc3/step_methods/metropolis.py index e0d049dba73..3daed1f7806 100644 --- a/pymc3/step_methods/metropolis.py +++ b/pymc3/step_methods/metropolis.py @@ -2,9 +2,10 @@ import numpy.random as nr import theano import scipy.linalg +import warnings from ..distributions import draw_values -from .arraystep import ArrayStepShared, ArrayStep, metrop_select, Competence +from .arraystep import ArrayStepShared, PopulationArrayStepShared, ArrayStep, metrop_select, Competence import pymc3 as pm from pymc3.theanof import floatX @@ -25,6 +26,11 @@ def __call__(self): return nr.normal(scale=self.s) +class UniformProposal(Proposal): + def __call__(self): + return nr.uniform(low=-self.s, high=self.s, size=len(self.s)) + + class CauchyProposal(Proposal): def __call__(self): return nr.standard_cauchy(size=np.size(self.s)) * self.s @@ -168,9 +174,7 @@ def astep(self, q0): @staticmethod def competence(var, has_grad): - if var.dtype in pm.discrete_types: - return Competence.COMPATIBLE - return Competence.INCOMPATIBLE + return Competence.COMPATIBLE def tune(scale, acc_rate): @@ -348,6 +352,7 @@ def competence(var): return Competence.IDEAL return Competence.INCOMPATIBLE + class CategoricalGibbsMetropolis(ArrayStep): """A Metropolis-within-Gibbs step method optimized for categorical variables. This step method works for Bernoulli variables as well, but it is not @@ -465,16 +470,135 @@ def competence(var): return Competence.COMPATIBLE return Competence.INCOMPATIBLE + +class DEMetropolis(PopulationArrayStepShared): + """ + Differential Evolution Metropolis sampling step. + + Parameters + ---------- + lamb : float + Lambda parameter of the DE proposal mechanism. Defaults to 2.38 / sqrt(2 * ndim) + vars : list + List of variables for sampler + S : standard deviation or covariance matrix + Some measure of variance to parameterize proposal distribution + proposal_dist : function + Function that returns zero-mean deviates when parameterized with + S (and n). Defaults to Uniform(-S,+S). + scaling : scalar or array + Initial scale factor for epsilon. Defaults to 0.001 + tune : bool + Flag for tuning the scaling. Defaults to True. + tune_interval : int + The frequency of tuning. Defaults to 100 iterations. + model : PyMC Model + Optional model for sampling step. Defaults to None (taken from context). + mode : string or `Mode` instance. + compilation mode passed to Theano functions + + References + ---------- + .. [Braak2006] Cajo C.F. ter Braak (2006). + A Markov Chain Monte Carlo version of the genetic algorithm + Differential Evolution: easy Bayesian computing for real parameter spaces. + Statistics and Computing + `link `__ + """ + name = 'DEMetropolis' + + default_blocked = True + generates_stats = True + stats_dtypes = [{ + 'accept': np.float64, + 'tune': np.bool, + }] + + def __init__(self, vars=None, S=None, proposal_dist=None, lamb=None, scaling=0.001, + tune=True, tune_interval=100, model=None, mode=None, **kwargs): + warnings.warn('Population based sampling methods such as DEMetropolis are experimental.' \ + ' Use carefully and be extra critical about their results!') + + model = pm.modelcontext(model) + + if vars is None: + vars = model.cont_vars + vars = pm.inputvars(vars) + + if S is None: + S = np.ones(sum(v.dsize for v in vars)) + + if proposal_dist is not None: + self.proposal_dist = proposal_dist(S) + else: + self.proposal_dist = UniformProposal(S) + + self.scaling = np.atleast_1d(scaling).astype('d') + if lamb is None: + lamb = 2.38 / np.sqrt(2 * S.size) + self.lamb = float(lamb) + self.tune = tune + self.tune_interval = tune_interval + self.steps_until_tune = tune_interval + self.accepted = 0 + + self.mode = mode + + shared = pm.make_shared_replacements(vars, model) + self.delta_logp = delta_logp(model.logpt, vars, shared) + super(DEMetropolis, self).__init__(vars, shared) + + def astep(self, q0): + if not self.steps_until_tune and self.tune: + # Tune scaling parameter + self.scaling = tune( + self.scaling, self.accepted / float(self.tune_interval)) + # Reset counter + self.steps_until_tune = self.tune_interval + self.accepted = 0 + + epsilon = self.proposal_dist() * self.scaling + + # differential evolution proposal + # select two other chains + ir1, ir2 = np.random.choice(self.other_chains, 2, replace=False) + r1 = self.bij.map(self.population[ir1]) + r2 = self.bij.map(self.population[ir2]) + # propose a jump + q = floatX(q0 + self.lamb * (r1 - r2) + epsilon) + + accept = self.delta_logp(q, q0) + q_new, accepted = metrop_select(accept, q, q0) + self.accepted += accepted + + self.steps_until_tune -= 1 + + stats = { + 'tune': self.tune, + 'accept': np.exp(accept), + } + + return q_new, [stats] + + @staticmethod + def competence(var, has_grad): + if var.dtype in pm.discrete_types: + return Competence.INCOMPATIBLE + return Competence.COMPATIBLE + + def sample_except(limit, excluded): candidate = nr.choice(limit - 1) if candidate >= excluded: candidate += 1 return candidate + def softmax(x): e_x = np.exp(x - np.max(x)) return e_x / np.sum(e_x, axis = 0) + def delta_logp(logp, vars, shared): [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared) diff --git a/pymc3/step_methods/slicer.py b/pymc3/step_methods/slicer.py index 3b5ec241cb0..f5f7d7e9c51 100644 --- a/pymc3/step_methods/slicer.py +++ b/pymc3/step_methods/slicer.py @@ -97,7 +97,7 @@ def astep(self, q0, logp): @staticmethod def competence(var, has_grad): if var.dtype in continuous_types: - if not var.shape: + if not var.shape or var.shape.ndim == 1: return Competence.PREFERRED return Competence.COMPATIBLE return Competence.INCOMPATIBLE diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py index 22eb7053ccd..8e12607d95d 100644 --- a/pymc3/tests/test_step.py +++ b/pymc3/tests/test_step.py @@ -9,7 +9,7 @@ from pymc3.step_methods import (NUTS, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, Metropolis, Slice, CompoundStep, NormalProposal, MultivariateNormalProposal, HamiltonianMC, - EllipticalSlice, smc) + EllipticalSlice, smc, DEMetropolis) from pymc3.theanof import floatX from pymc3 import SamplingError from pymc3.distributions import ( @@ -318,7 +318,7 @@ def test_mv_proposal(self): class TestCompoundStep(object): - samplers = (Metropolis, Slice, HamiltonianMC, NUTS) + samplers = (Metropolis, Slice, HamiltonianMC, NUTS, DEMetropolis) @pytest.mark.skipif(theano.config.floatX == "float32", reason="Test fails on 32 bit due to linalg issues") @@ -393,6 +393,22 @@ def kill_grad(x): assert isinstance(steps, Slice) +class TestPopulationSamplers(object): + def test_checks_population_size(self): + """Test that population samplers check the population size.""" + steppers = [ + DEMetropolis + ] + with Model() as model: + n = Normal('n', mu=0, sd=1) + for stepper in steppers: + step = stepper() + with pytest.raises(ValueError): + trace = sample(draws=100, chains=1, step=step) + trace = sample(draws=100, chains=4, step=step) + pass + + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") class TestNutsCheckTrace(object): def test_multiple_samplers(self):