Gaussian prior implementation

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!)
 ===========================================================

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

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'

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)