diff --git a/README.rst b/README.rst index 5dd9abb5..33877149 100644 --- a/README.rst +++ b/README.rst @@ -1,3 +1,22 @@ +Gaussian priors for MultiNest sampling: +------------------- + +The standard method of putting a Gaussian prior on sampling parameters (i.e. adding a Gaussian likelihood) leads to an incorrect +estimate of the evidence with MultiNest. This branch contains an implementation of a Gaussian prior in ``montepython/prior.py`` +(adapted from https://github.com/JohannesBuchner/MultiNest/blob/master/src/priors.f90). + +The mean and standard deviation of the Gaussian are defined in the ``.param`` file. For example: + + data.parameters['h'] = [0.68, 0.6, 0.8, 0.2, 1, 'cosmo', 'gaussian', 0.68, 0.001] + +puts a Gaussian prior on h with mu=0.68 and sigma=0.001. If we use + + data.parameters['omega_cdm'] = [0.1120, -1,-1, 0.0016,1, 'cosmo'] + +we assume a flat prior in the given parameter range. + +Note: At the moment the Gaussian prior only works with the MultiNest sampler, but it should be straightforward to implement in PolyChord as well. + =========================================================== Monte Python, a Monte Carlo Markov Chain code (with Class!) =========================================================== diff --git a/montepython/MultiNest.py b/montepython/MultiNest.py index c67d3e95..d224b779 100644 --- a/montepython/MultiNest.py +++ b/montepython/MultiNest.py @@ -115,10 +115,10 @@ def initialise(cosmo, data, command_line): # Check that all the priors are flat and that all the parameters are bound is_flat, is_bound = sampler.check_flat_bound_priors( data.mcmc_parameters, varying_param_names) - if not is_flat: - raise io_mp.ConfigurationError( - 'Nested Sampling with MultiNest is only possible with flat ' + - 'priors. Sorry!') + # if not is_flat: + # raise io_mp.ConfigurationError( + # 'Nested Sampling with MultiNest is only possible with flat ' + + # 'priors. Sorry!') if not is_bound: raise io_mp.ConfigurationError( 'Nested Sampling with MultiNest is only possible for bound ' + @@ -169,7 +169,7 @@ def initialise(cosmo, data, command_line): if not param in NS_param_names: NS_param_names.append(param) data.NS_param_names = NS_param_names - + # Caveat: multi-modal sampling OFF by default; if requested, INS disabled try: if data.NS_arguments['multimodal']: @@ -273,7 +273,7 @@ def loglike(cube, ndim, *args): for i, name in enumerate(derived_param_names): cube[ndim+i] = data.mcmc_parameters[name]['current'] return lkl - + #FK: recover name of base folder and remove entry from dict before passing it # on to MN: base_dir = data.NS_arguments['base_dir'] @@ -356,10 +356,10 @@ def from_NS_output_to_chains(folder): if line.strip()[0] == '#': continue - # These lines allow MultiNest to deal with fixed nuisance parameters + # These lines allow MultiNest to deal with fixed nuisance parameters sigma = float(line.split(',')[3].strip()) if sigma == 0.0: - #If derived parameter, keep it, else discard it: + #If derived parameter, keep it, else discard it: paramtype = line.split(',')[5].strip()[1:-2] if paramtype != 'derived': continue diff --git a/montepython/data.py b/montepython/data.py index 1d6fd542..eecb9889 100644 --- a/montepython/data.py +++ b/montepython/data.py @@ -169,7 +169,7 @@ def __init__(self, command_line, path): :rtype: dict """ - # Arguments for PyPolyChord + # Arguments for PyPolyChord self.PC_param_names = [] self.PC_arguments = {} """ @@ -1119,7 +1119,7 @@ def __init__(self, array, key): self['initial'] = array[0:4] self['scale'] = array[4] - self['role'] = array[-1] + self['role'] = array[5] self['tex_name'] = io_mp.get_tex_name(key) if array[3] == 0: self['status'] = 'fixed' diff --git a/montepython/prior.py b/montepython/prior.py index 24171579..9e0c55ad 100644 --- a/montepython/prior.py +++ b/montepython/prior.py @@ -8,6 +8,8 @@ import random as rd from copy import deepcopy import io_mp +from scipy.special import erfcinv +import numpy as np class Prior(object): """ @@ -79,7 +81,7 @@ def draw_from_prior(self): within_bounds = calue_within_prior_range(value) return value - + def value_within_prior_range(self, value): """ Check for a value being in or outside the prior range @@ -112,5 +114,13 @@ def map_from_unit_interval(self, value): which should have been previously checked with :func:`is_bound` """ - return (self.prior_range[0] + + if self.prior_type == 'flat': + return (self.prior_range[0] + value * (self.prior_range[1] - self.prior_range[0])) + # Gaussian prior adapted from https://github.com/JohannesBuchner/MultiNest/blob/master/src/priors.f90 + if self.prior_type == 'gaussian': + if (value <= 1e-16) or ((1 - value) <= 1e-16): + gaussian_prior = -1e-32 + else: + gaussian_prior = self.mu + self.sigma * np.sqrt(2) * erfcinv(2*(1-value)) + return(gaussian_prior)