Parallel-Tempered MCMC (ligo.skymap.extern.ptemcee)

A parallel-tempered version of emcee, adapted from https://github.com/willvousden/ptemcee. See ptemcee License.

class ligo.skymap.extern.ptemcee.InterruptiblePool(processes=None, initializer=None, initargs=(), **kwargs)

A modified version of multiprocessing.pool.Pool that has better behavior with regard to KeyboardInterrupts in the map() method.

Parameters:

• processes – (optional) The number of worker processes to use; defaults to the number of CPUs.

• initializer – (optional) Either None, or a callable that will be invoked by each worker process when it starts.

• initargs – (optional) Arguments for initializer; it will be called as initializer(*initargs).

• kwargs – (optional) Extra arguments. Python 2.7 supports a maxtasksperchild parameter.

map(func, iterable, chunksize=None)

Equivalent of map() built-in, without swallowing KeyboardInterrupt.

Parameters:

• func – The function to apply to the items.

• iterable – An iterable of items that will have func applied to them.

class ligo.skymap.extern.ptemcee.MPIPool(comm=None, debug=False, loadbalance=False)

A pool that distributes tasks over a set of MPI processes. MPI is an API for distributed memory parallelism. This pool will let you run emcee without shared memory, letting you use much larger machines with emcee.

The pool only support the map() method at the moment because this is the only functionality that emcee needs. That being said, this pool is fairly general and it could be used for other purposes.

Contributed by Joe Zuntz.

Parameters:

• comm – (optional) The mpi4py communicator.

• debug – (optional) If True, print out a lot of status updates at each step.

• loadbalance – (optional) if True and ntask > Ncpus, tries to loadbalance by sending out one task to each cpu first and then sending out the rest as the cpus get done.

bcast(*args, **kwargs)

Equivalent to mpi4py bcast() collective operation.

close()

Just send a message off to all the pool members which contains the special _close_pool_message sentinel.

is_master()

Is the current process the master?

map(function, tasks)

Like the built-in map() function, apply a function to all of the values in a list and return the list of results.

Parameters:

• function – The function to apply to the list.

• tasks – The list of elements.

wait()

If this isn’t the master process, wait for instructions.

class ligo.skymap.extern.ptemcee.Sampler(nwalkers, dim, logl, logp, ntemps=None, Tmax=None, betas=None, threads=1, pool=None, a=2.0, loglargs=[], logpargs=[], loglkwargs={}, logpkwargs={}, adaptation_lag=10000, adaptation_time=100, random=None, vectorize=False)

A parallel-tempered ensemble sampler, using EnsembleSampler for sampling within each parallel chain.

Parameters:

• nwalkers – The number of ensemble walkers at each temperature.

• dim – The dimension of parameter space.

• betas – (optional) Array giving the inverse temperatures, \(\beta=1/T\), used in the ladder. The default is chosen according to default_beta_ladder() using ntemps and Tmax.

• ntemps – (optional) If set, the number of temperatures to use.

• Tmax – (optional) If set, the maximum temperature for the ladder.

• logl – The log-likelihood function.

• logp – The log-prior function.

• threads – (optional) The number of parallel threads to use in sampling.

• pool – (optional) Alternative to threads. Any object that implements a map method compatible with the built-in map will do here. For example, multi.Pool will do.

• a – (optional) Proposal scale factor.

• loglargs – (optional) Positional arguments for the log-likelihood function.

• logpargs – (optional) Positional arguments for the log-prior function.

• loglkwargs – (optional) Keyword arguments for the log-likelihood function.

• logpkwargs – (optional) Keyword arguments for the log-prior function.

• adaptation_lag – (optional) Time lag for temperature dynamics decay. Default: 10000.

• adaptation_time – (optional) Time-scale for temperature dynamics. Default: 100.

• vectorize – (optional) If True, logl and logp are expected to accept a list of position vectors instead of just one. Note that pool will be ignored if this is True. (default: False)

property acceptance_fraction

Matrix of shape (Ntemps, Nwalkers) detailing the acceptance fraction for each walker.

property acor

Returns a matrix of autocorrelation lengths for each parameter in each temperature of shape (Ntemps, Ndim).

property beta_history

Matrix of inverse temperatures; shape (Ntemps, Nsteps).

property betas

Returns the current inverse temperature ladder of the sampler.

property chain

Returns the stored chain of samples; shape (Ntemps, Nwalkers, Nsteps, Ndim).

property flatchain

Returns the stored chain, but flattened along the walker axis, so of shape (Ntemps, Nwalkers*Nsteps, Ndim).

get_autocorr_time(window=50)

Returns a matrix of autocorrelation lengths for each parameter in each temperature of shape (Ntemps, Ndim).

Parameters:

window – (optional) The size of the windowing function. This is equivalent to the maximum number of lags to use. (default: 50)

log_evidence_estimate(logls=None, fburnin=0.1)

Thermodynamic integration estimate of the evidence for the sampler.

Parameters:

• logls – (optional) The log-likelihoods to use for computing the thermodynamic evidence. If None (the default), use the stored log-likelihoods in the sampler. Should be of shape (Ntemps, Nwalkers, Nsamples).

• fburnin – (optional) The fraction of the chain to discard as burnin samples; only the final 1-fburnin fraction of the samples will be used to compute the evidence; the default is fburnin = 0.1.

Return (logZ, dlogZ)(logZ, dlogZ):

Returns an estimate of the log-evidence and the error associated with the finite number of temperatures at which the posterior has been sampled.

For details, see thermodynamic_integration_log_evidence.

property loglikelihood

Matrix of log-likelihood values; shape (Ntemps, Nwalkers, Nsteps).

property logprobability

Matrix of logprobability values; shape (Ntemps, Nwalkers, Nsteps).

property ntemps

The number of temperature chains.

property random

Returns the random number generator for the sampler.

reset(random=None, betas=None, time=None)

Clear the time, chain, logposterior, loglikelihood, acceptance_fraction, tswap_acceptance_fraction stored properties.

run_mcmc(*args, **kwargs)

Identical to sample, but returns the final ensemble and discards intermediate ensembles.

sample(p0=None, iterations=1, thin=1, storechain=True, adapt=False, swap_ratios=False)

Advance the chains iterations steps as a generator.

Parameters:

• p0 – The initial positions of the walkers. Shape should be (ntemps, nwalkers, dim). Can be omitted if resuming from a previous run.

• iterations – (optional) The number of iterations to perform.

• thin – (optional) The number of iterations to perform between saving the state to the internal chain.

• storechain – (optional) If True store the iterations in the chain property.

• adapt – (optional) If True, the temperature ladder is dynamically adapted as the sampler runs to achieve uniform swap acceptance ratios between adjacent chains. See arXiv:1501.05823 for details.

• swap_ratios – (optional) If True, also yield the instantaneous inter-chain acceptance ratios in the fourth element of the yielded tuple.

At each iteration, this generator yields

• p, the current position of the walkers.

• logpost, the current posterior values for the walkers.

• loglike, the current likelihood values for the walkers.

• ratios, the instantaneous inter-chain acceptance ratios (if requested).

property time

Returns the current time, in iterations, of the sampler.

property tswap_acceptance_fraction

Returns an array of accepted temperature swap fractions for each temperature; shape (ntemps, ).

ligo.skymap.extern.ptemcee.default_beta_ladder(ndim, ntemps=None, Tmax=None)

Returns a ladder of \(eta \equiv 1/T\) under a geometric spacing that is determined by the arguments ntemps and Tmax. The temperature selection algorithm works as follows:

Ideally, Tmax should be specified such that the tempered posterior looks like the prior at this temperature. If using adaptive parallel tempering, per arXiv:1501.05823, choosing Tmax = inf is a safe bet, so long as ntemps is also specified.

Parameters:

• ndim – The number of dimensions in the parameter space.

• ntemps – (optional) If set, the number of temperatures to generate.

• Tmax – (optional) If set, the maximum temperature for the ladder.

Temperatures are chosen according to the following algorithm:

• If neither ntemps nor Tmax is specified, raise an exception (insufficient information).

• If ntemps is specified but not Tmax, return a ladder spaced so that a Gaussian posterior would have a 25% temperature swap acceptance ratio.

• If Tmax is specified but not ntemps:

  • If Tmax = inf, raise an exception (insufficient information).

  • Else, space chains geometrically as above (for 25% acceptance) until Tmax is reached.

• If Tmax and ntemps are specified:

  • If Tmax = inf, place one chain at inf and ntemps-1 in a 25% geometric spacing.

  • Else, use the unique geometric spacing defined by ntemps and Tmax.