geometric nested sampling

The README below as well as the package API can be found at: gns read the docs.

Version

Version: 1.01

Overview

Nested sampling algorithms are a family of algorithms designed to sample from posterior probability distributions associated with e.g. Bayesian inference. The overall idea in nested sampling is to convert an n-dimensional integral over the parameter space into a single dimensional integral over 'prior volume' and equivalently sample from this domain. The purpose of the geometric nested sampling algorithm presented here is to efficiently sample the distributions of parameters which exhibit certain geometrical properties e.g. circular, toroidal or spherical. Such parameters are rife in the astrophysics community, where many time- or angle-periodic parameters are present.

The algorithm

Can run Markov Chain Monte Carlo (MCMC) Metropolis Hastings (MH) nested sampling algorithm in 'vanilla' (see nested sampling and ns with MCMC) and/or 'geometric' (gns paper) sampling modes. See the appropriate papers for detailed explanations of how these algorithms operate.

In vanilla mode, the trial points are sampled in the physical (native) space of the input parameters. In geometric mode (main point of the algorithm), parameter values represent points on a geometric object (e.g. circle, torus or sphere), and trial points are picked by sampling from spaces parameterised by Cartesian coordinates. If the sampled points do not lay on the original object (in general they will not), then they are projected back onto the object before being transformed back to the native space for evaluation of the MH acceptance ratio.

The sampling mode for each parameter is specified in paramGeomList passed through setupDict to nestedRun() ("vanilla", "wrapped" or "sphere"). "wrapped" parameters use a wrapped trial distribution. "sphere" parameters are transformed into their Euclidean representation and sampled in this space. Geometric sampling works best for parameters with 'periodic' sampling domains e.g. angles, as sampling in geometric space results in wrapping parameter values at end of domains. The sampler also works well for parameters which exhibit 'spherical' properties. Note that when sampling spherical parameters, two physical parameters are required for each geometric transformation, i.e. the spherical parameters must come in pairs. See tests.py for examples of how to specify what parameters are sampled in what way and the below section showing paramGeomList examples.

Main output of geometric nested sampling run is file of sampled points and their weights. These are stored in text_output/file_name.txt, whose format follows that explained in getdist RTD. Similarly, metadata files are also generated for each run e.g. file_name.paramnames, see getdist RTD for more details on these files. The format of the text outputs is such that they can be plotted very easily using getdist, a kernel density esimate plotter.

Toy models

To run with toy models (defined and explained in toy_models.py), go to tests/ and run tests.py, which calls the function runTests(), which can be found in the tests.py file. Plots for comparison can be found in image_output/toy_models_empirical_getdist/. Text output can be found in text_output/. The default toy model used when main.py is ran as is when cloned from this repo, should be compared with files image_output/toy_models_empirical_getdist/circle.png (in particular, the MHWG curve on this plot). Files created from running main.py for this toy example are saved as th_Sc_mhwg* files in text_output/ and image_output/toy_models_empirical_getdist. Further toy models can be run by adding them to tests.py and adjusting the necessary parameters (see the commented out example for how to configure the parameters). The following plot was produced from running the default circle toy model, and shows the kernel density esimate of samples obtained by the geometric nested sampler from the von Mises distribution.

Examples of setupDict objects used to specify parameters of a nested sampling run can be found in setup_dicts.py. 'verbose' corresponds to whether to output text during the sampling run, 'trapezoidalFlag' specifies whether to use the trapezoidal rectangular quadrature techniques for approximating the integral, 'ZLiveType' specifies how to determine the final contribution to the nested sample run, and shouldn't be significant if the algorithm is ran for a decent amount of time, 'terminationType' is a similar parameter, as is 'terminationFactor', 'sampler' is the type of sampler used i.e. Blind nested sampling, which involves randomly sampling values from the unit hypercube (B) Metropolis Hastings nested sampling (MH) geometric Metropolis Hastings nested sampling (MHWG), or use MultiNest if it is installed (MN). 'outputFile': is the location to store the text output from the runs, 'space' corresponds to whether the likelihood or log likelihood is used in calculations (and thus the corresponding function should be passed to the nested sampling algorithm through the appropriate arguments, see nestedRun()). 'paramGeomList': is a list of size of the dimensionality of the problem, and specifies how each parameter should be treated by the geometric nested sampler: "vanilla", "wrapped" or "sphere" (must also be specified for the MH algorithm, which only accepts 'vanilla' values) Note for the toy models, none of these really need to be changed from their default in tests.py apart from paramGeomList, only the models chosen and the prefixes/suffixes used to index them need to be altered.

The toy models (see https://arxiv.org/abs/1905.09110) are designed to make use of the key features of geometric nested sampling, i.e. to take advantage of the way it treats parameters defined on curved manifolds. The three toy models which feature in tests.py by default are a circular distribution (von Mises distribution), a toroidal distribution (two-dimensional von Mises distribution), and 6 spheres with Kent distributions defined on them (for more on this see gns paper and Kent paper). Many more are available in toy_models.py, all that needs to be configured is there paramGeomList values when added to tests.py, and the relevant suffixes/prefixes to index them.

Example

The following script produces samples from the von Mises distribution shown in the above image, and plots these samples using getdist. Note to produce the plots at the end of the example, getdist must be installed. This will be done automatically if gns was pip installed (see Installation section below) with the plotting requirements specified, i.e. pip install GNS[plotting]==1.0. Otherwise getdist can be installed with the command pip install getdist.

import numpy as np
import os
import gns.nested_run as nested_run
import gns.plotting as plotting
import gns.toy_models as toy_models

#see previous section for what different pairs in setupDict represent. 
setupDict = {'verbose':True, 'trapezoidalFlag': False, 'ZLiveType':'average X', 'terminationType':'evidence', 'terminationFactor': 0.1, 'sampler':None, 'outputFile':None, 'space':'log'}

shape = 'circle' #run circle toy model, simple von Mises distribution
sampler = 'MH WG' #geometric nested sampler

#make sure directories to store output exist, and if not, create them
textDir = './text_output/'
plotDir = './image_output/'
if not os.path.exists(textDir):
  os.mkdir(textDir)
if not os.path.exists(plotDir):
  os.mkdir(plotDir)

outputFile = textDir + 'example' #for text output
plotFile = plotDir + 'example' #for image output

#get prior and likelihood functions from toy models for circle toy model
paramNames, priorParams, LhoodParams = toy_models.getToyHypersGeom(shape)
priorFunc, logPriorFunc, invPriorFunc, _, LLhoodFunc = toy_models.getToyFuncs(priorParams, LhoodParams)

#can manually set paramNames e.g.
paramNames = [r"\phi"] #label plot

nDims = 1 #number of parameters in inference
#array with first element giving lower bound of prior support, second giving upper bound, and third giving upper - lower. shape (3, nDims)
targetSupport = np.array([0., 2. * np.pi, 2. * np.pi]).reshape(3,1) 

#See next section for examples of paramGeomList. For simple one-dimensional circular example considered here, just need to consider one 'wrapped' parameter
setupDict['paramGeomList'] = ['wrapped']

#output files
setupDict['outputFile'] = outputFile
setupDict['sampler'] = sampler

#run sampler
nested_run.NestedRun(priorFunc, invPriorFunc, None, paramNames, targetSupport, setupDict, LLhoodFunc)

#plot data using getdist
#NOTE, calls to output.writeParamNames() and output.writeRanges() are made in nested_run to create .paramnames and .ranges files required by getdist
plotting.callGetDist([outputFile], plotFile + '.png', nDims, ['mhwg'])

Running the algorithm with custom likelihood/prior functions

If custom priors/ likelihoods are to be used, one needs to call the nestedRun() function in the gns.nested_run module, and pass the required likelihood and prior functions. All of these functions must take as an argument a numpy array with shape (nLive, nDims) where nLive is the number of live points in the nested run, and is defined in the functions NestedRunLinear and NestedRunLog which are in the nested_run.py file. nDims is the dimension