A wrapper for use of SSPR in nested sampling packages such as PolyChord and Multinest
Stochastic superpositional posterior repartitioning used in conjunction with Bayesian inference engines such as PolyChord and MultiNest.
What is Bayesian inference
Bayesian inference is a suitable framework for model fitting. You give it a theory and some data and it tells you how well the theory fits the data, while also telling you what the theory's parameters should be for the fit to be the best.
The preferred way is to use
pip install super-nest
though other packages will be developed as needed.
If you know what you're doing you can clone this repository and install manually, using
cd super-nest && python3 setup.py
You also need to have either PolyChord or MultiNest (or both) installed.
This is planned as a simple convenient package that you use with a sampler.
Suppose that you had a set-up with
PolyChord that involved a prior
prior a log-likelihood function,
loglike all of which
operated on an
nDims dimensional space.
Often you'd be using a uniform prior. You're sure that using a different prior would cause an imprint, which is fine if that prior was based on a previous run, but not so if you just used it as a hunch. As hunch is used liberally here, you could have an intuition, or you could have done a run using a different model and extrapolated the posterior.
Previously all such information could not be used. Nested sampling would return all of that information to you and you wouldn't be able to tell, if this is actually what the data suggests, or that you simply gave it too much info out of thin air.
Well, now you can use that information and get your sampling to run faster, but to produce the output you would have gotten had you used a uniform prior (and a lot more live points).
Suppose you had a unfirom prior
pi and likelihood
l. To use them
in nested sampling, you need a prior quantile, and the logarithm of
the likelihood. You need to pass both to the sampler, e.g.
# This is pseudo-code. settings = PolyChordSettings(nDims, nDerived, ...) run_polychord(log(l), functional_inverse(cumulative_dist(pi)), nDims, nDerived ...)
if the prior is uniform (a function that maps to a
constant), you expect the quantile to be dependent on the boundaries
of your prior space. If you have a box with corners at
then the uniform prior quantile is:
def quantile(cube): return a + (b-a) * cube
Let's say you have a hunch that the posterior would be a Gaussian at
mu with covariance matrix
Sigma. How do you formulate it in such a
way that supernest can understand?
Well, for starters, the prior quantile would look like
def gaussian_quantile(cube): def _erf_term(d, g): return erf(d * sqrt(1 / 2) / g) sigma = diag(m.cov) da = _erf_term(a - mu, sigma) db = _erf_term(b - mu, sigma) ret = erfinv((1 - cube) * da + cube * db) return mu + sqrt(2 / beta) * sigma * ret
which is already quite a mouthful. If you tried to use this quantile directly, you would have a few problems:
- Your posterior will have a Gaussian imprint at
Sigma, whether or not you want it.
- If most of your posterior is far away, it would have been picked up by the uniform prior, but not the one you have.
- You would get the wrong (larger) evidence than you would with a uniform prior.
So what do you do?
Creating a proposal.
You create a consistent proposal. This means that you need to change the likelihood specifically, you need to make it so that the product of the prior pdf and the new likelihood is the same as the original uniform prior times the original
So what do you do for a Gaussian? Well, it's easy:
def logl_gaussian(theta): ll=0 def ln_z(t, b): ret = - b * (t - mu) ** 2 / 2 / sigma ** 2 ret -= log(pi * sigma ** 2 / 2 / b) / 2 db = _erf_term(b - mu, sigma) da = _erf_term(a - mu, sigma) ret -= log(db - da) return ret ll -= log(b-a).sum() ll -= ln_z(theta, beta).sum() return ll
And then you pass it and it works? Wrong. This will help you with the wrong evidence, but you would still get an imprint in the posterior, and sampling far away would still be an issue.
Using this package.
To solve the problem you need to put both your original prior and likelihood and the newly created gaussian proposa. into a superposition. This what this package does, and this is what this package is.
So to use the things that we've defined previously, we should do the following: (I'm assuming polyChord, but this works with any nested sampler).
from polychord import run_polychord from polychord.settings import PolyChordSettings from supernest import superimpose def uniform_quantile(cube): ... def uniform_like(theta): ... return log(l), derived def gaussian_quantile(cube): ... # In later versions, you will have a function that generates the correctsion for you. def gaussian_like(theta): ... # Note taht you can have as many proposals as you like. n_dims, prior, log_like = superimpose([(uniform_prior, uniform_like), (gaussian_prior, gaussian_like), ...], nDims) ... # Be wary of using grade_dims and grade_frac. If you have n models, you should # grade_dims.append(n) # grade_frac.append([1.0]*n) # before passing them to settings. settings = PolyChordSettings(n_dims,nDerived, ...) ... run_polychord(prior, n_dims, nDerived, log_like, ...)
What do I get in return for my work?
A run-time reduction starting at 37x. A precision increase of about 100x. The posterior you get should contain all the usual stuff, plus some extra parameters at the end, which tell you how good your proposals were. These should be interpreted as probabilities that the proposal described the posterior.
Is this the best that superpositional mixtures have to offer.
Not at all. We're writing a purpose built sampler that eliminates all outward appearances of using superposition under the hood, while also making use of some clever tricks, that you normally can't do. Effectively we're creating a quantum computer of nested samplers.
Is there a neater way?
Yes. If you know exactly what you want, you can (optionally) use the
super_nest.framework module. This is a more class oriented interface that's much nicer to deal with, as it infers the dimensionality of the problem and is far simpler to use than the non-OOP direct interface.
So in that case, you would do something like
from super_nest.framework.gaussian_models.uniform import BoxUniformModel mdl = BoxUniformModel((a, b), mu, cov, file_root='Uniform')
which creates a model with a box uniform prior which has corners at
b, and a single Gaussian peak at
mu with covariance
Why would you want to use this? Well, how about if you wanted to create a mixture of this model, along with a custom model, that you defined as a class?
from super_nest.framework.polychord_model import Model class MyCustomModel(Model): def __init__(self, *args, **kwargs): # get some stuff you need pass def log_likelihood(self, theta): return myCustomLogLike(theta) def quantiel(self, hypercube): return myPriorQuantile(hypercube) @property def dimensionality(self): return nDims @property def num_derived(self): return nDerived
Which is admittedly more verbose, than just using base
super_nest.superimpose. The reason why you'd want to do that, is
that this will automatically track the dimensionality for you, create
loglike and prior quantile that correspond to what you
want, and allows you to do some non-standard posterior repartitioning,
that would otherwise be a headache to deal with.
TL;DR, you could just do things like
from super_nest.framework.mixture import StochasticMixtureModel from super_nest.framework.gaussian_models import (GaussianPeakedPrior, Uniform, PowerPosteriorPrior) uniform = Uniform((a,b), mu, cov, ...) gaussian = GaussianPeakedPrior((a,b), mu, cov, ...) ppr = PowerPosteriorPrior((a,b), mu, cov,...) mix = StochasticMixtureModel([uniform, gaussian, ppr]) ... mix.nested_sample()
which for larger projects where the off by one error doesn't always cause a segfault, but is something that you discover several HPC hours later, is a necessity. I created this, despite having a decent-enough interface to PolyChord, because the situation of a Gaussian proposal mixed with some arbitrary code is a common use case. And in some circumstances (e.g. where you have no idea what the confidence intervals might be), you could prefer to use power posterior repartitioning, inside of a stochastic mixture.
You can do it with just
super_nest.superimpose, but it would be a
much bigger hassle. The domain of Bayesian inference has many common
patters, so it really does pay to use a slightly more complicated
system, but not to worry about "did I get the linear transformation
right?" kinds of questions. This is MIT licensed, so why bother
implementing something in your own code, when you can just create a PR
(I promise I don't bite).
Unfortunately this framework only works with PolyChord (yet), but work
is in place to make it universal. The idea is that different samplers
have different use-cases, and ideally the
choose the right tool for the right job. Specifically, there's work
planned to include support for
PyMultiNest, as it's much faster for
low-dimensional inference problems (it scales exponentially, so low
really means no more than five parameters). In some cases, when the
log-likelihood evaluation is very slow, it would make sense to use
dynesty, while pure-python has its drawbacks (the
only thing this language is good at is racking massive technical debt,
that you hope never to reap),
dynesty is a dynamic sampler,
meaning that it can potentially outperform any of the ones on the list
dypolychord, which is a very well-made package).
It will likely be the chief influence on the API for the next generation dynamic superpositional trans-dimensional sampler, so getting to know this API might be worth the extra time.
Does this support the X sampler.
Yes. This is a very thin module, and while I have assumed a
PolyChord-like calling convention (also
dynesty), you can use it
with other samplers that support using python functions as callbacks.
Very similar to PolyChord, (and also quite popular). Use as follows:
from super_nest import superimpose n_dims, my_prior_quantile, my_likelihood = superimpose( [(prior1, like1), (proposal_prior1, proposal_like1)], nDims) solve(LogLikelihood=my_likelihood, Prior=my_prior_quantile, n_dims=n_dims)
Has the same standard interface.
IF you want to use
Multinest, you should wait for a FORTRAN version of
this package to use directly with MultiNest. If you want to use
pymultinest a hard
requirement for now.
What is this useful for
Suppose that you were running a complex Bayesian inference, e.g. Cobaya. You have a choice of sampler, e.g. Monte-Carlo, Metropolis Hastings or you could choose to use Nested Sampling, among a family of other inference methods. If you choose the former, you get a good idea of what the model parameters should be quick, but you have no idea how good the fit is, because MCMC and MH don't evaluate the evidence. You think to use Nested Sampling, but that takes too long, and you can't give it more information to run faster... At least not without this package.
This is a thin wrapper around both PolyChord and Multinest's code, that allows you to specify a proposal distribution. For now it's mainly a multivariate correlated Gaussian, but other distributions are being planned.
How much faster/more precise can this make inference?
Our preliminary academic benchmarks showed a runtime reduction by a factor of 20. Your mileage may vary, but if you're willing to sacrifice some precision, you can make it go even faster.
If you really want the extra precision, you can expect an uplift by two orders of magnitude, if you started with a uniform prior (as we often do in nested sampling) and used the posterior chains to generate the distribution.
License - LGPLv3
It's a permissive license that's GPL compatible. if you want a deeper understanding: (https://www.gnu.org/licenses/lgpl-3.0.en.html)[go here].
This is a python package. You don't normally need much beyond creating pull requests, for features that you think should be in.
I don't care about Pep8. It's a misguided collection of near-sighted practices, that are emblematic of the problems python has. Namely, that it's an interpreted language that's very hard to read for a computer, making python slower than it has to be.
So this means that
- if you have a proposal, and it's not pep8. It's fine.
- if you have a proposal and it's pep8, it's also fine.
- if you have a proposal that's just about making the code pep8. It's fine.
I don't bite.
Release history Release notifications | RSS feed
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Hashes for super_nest-1.0.0-py3-none-any.whl