baudren / montepython_public

Public repository for the Monte Python Code
MIT License
65 stars 114 forks source link

How to use Gaussian priors with MultiNest? #112

Closed fkoehlin closed 6 years ago

fkoehlin commented 6 years ago

Dear all,

I understand that once you use MultiNest as a sampler for MontePython all priors on cosmology and nuisance parameters are in general interpreted as flat priors within the given finite ranges (i.e. their settings for the mean and variance are ignored). However, is it possible to still enforce Gaussian priors on 'cosmo' and/or 'nuisance' parameters by using settings like

a) data.parameters['H0'] = [72., 60., 80., 1.0, 1, 'cosmo', 'gaussian', 73., 2.00] (being aware of the potential bug as pointed out in issue #55) when using MultiNest (and could I just replace 'cosmo' with 'nuisance' then, too)?

In addition to a) I've found in the documentation for MP (version 2.2.0; page 25) a section on Gaussian priors for nuisance parameters which need to be defined within the 'likelihood_name.data' file like that:

b) 'likelihood_name.nuisanceparamprior_center = X' and 'likelihood_name.nuisanceparamprior_variance = Y'

So my questions in short:

1) Is there a way to enforce Gaussian priors when using MultiNest? 2) What's the difference between a) and b) (if a) would work with 'nuisance', too)?

Many thanks for your answers in advance!

Cheers,

Fabian

brinckmann commented 6 years ago

Hi Fabian,

I do not believe MultiNest can handle priors in the first way you suggest, as the only information relevant to MultiNest are your parameter ranges and MultiNest settings.

However, you can add a prior as you suggest in method b, e.g. by following the example of the "hst" likelihood.

To specifically address your questions: 1) modify 'hst' to use the center and variance you desire for H0, or create a new likelihood modifying all instances of hst to new_likelihood, i.e. foldername, .data filename and contents and the class definition in the init file. 2) The former [a)] changes the sampling distribution. Instead of steps being generated from a flat distribution, they will be generated from a Gaussian distribution. The latter [b)] accepts/rejects steps according to the new Gaussian prior likelihood, just as if it was any other likelihood. It would work with a nuisance parameter, you just need to correctly call your nuisance parameter in the init file, i.e. replace h in the hst likelihood with nuisance_parameter = data.mcmc_parameters['nuisance_parameter']['current']*data.mcmc_parameters['nuisance_parameter']['scale']

Best, Thejs

fkoehlin commented 6 years ago

Hi Thejs,

Thanks for your fast reply! I fully get your explanation for [a)] but I'm confused now with the "hst" likelihood example and 'my' method [b)]: do you expect the recipe given in documentation for MP (version 2.2.0; page 25), i.e. for example defining two nuisance parameters like

likelihood.use_nuisance = ['N1', 'N2'] likelihood.N1_prior_center = 1. likelihood.N1_prior_variance = 2. likelihood.N2_prior_center = 1. likelihood.N2_prior_variance = 2.

directly in the 'likelihood.data'-file of the main likelihood for several/a nuisance parameter/s to work with MultiNest? Or do I further need to modify sth. in the main likelihood (e.g. 'chi2 += (N2-center)**2/variance')?

And how does your "hst" example fit in here? In the sense of having to define a 'nuisance prior-likelihood' for all nuisances and then sampling over those explicitly by adding them to the list of likelihoods in the parameter-file? Sorry for the confusion and thanks in advance for your help!

Cheers, Fabian

brinckmann commented 6 years ago

The example you mention from the documentation is just a variant of a simple likelihood, similar to the hst likelihood (except it is a prior on two nuisance parameters rather than a cosmological parameter). I just used hst as an example because it is a typical example of a likelihood working as a Gaussian prior on some parameter.

You can proceed in exactly the same way, either by creating a new likelihood or by modifying your main likelihood. In either case you would need to add the contribution to chi2 as you say.

Best, Thejs

fkoehlin commented 6 years ago

Hi Thejs,

Thanks again for your help. My confusion was based on reading over that the example given in the doc was for a specific likelihood ('sn' in that case). Hence, I thought that every likelihood would inherit the '..._prior_center' and '..._prior_variance' keywords. Long story short: I've implemented now my own additions to the chi^2 for the nuisance parameters used in my likelihood and that allows me to define Gaussian priors on them, even when using MultiNest.

In case this issue might pop up in future searches, here's what I did:

Add to your likelihood.data file the following lines:

# if you want to enforce Gaussian priors on some/all NUISANCE parameters, set flag to True
# if set to False lists below are ignored!
likelihood.use_gaussian_prior_for_nuisance = True
# add here all NUISANCE parameters for which you would like to define Gaussian priors:
# name must match exactly to the nuisances defined above and must appear in list 'use_nuisance'!
likelihood.gaussian_prior_name = ['nuisance_1', 'nuisance_2', ..., 'nuisance_N']
# supply here the central values of the Gaussians (keep the order! no double checks!!!)
likelihood.gaussian_prior_center = [1.0, 1.1, ..., 1.2]
# supply here the std wrt. the center (again keep the order!)
likelihood.gaussian_prior_sigma = [0.1, 0.2, ..., 0.5]

Add to your likelihood (__init__.py) within the 'loglkl'-function just before the return:

chi2 = calc_chi2(data, model, covariance)

# enforce Gaussian priors on NUISANCE parameters if requested:
if self.use_gaussian_prior_for_nuisance:
    for idx_nuisance, nuisance_name in enumerate(self.gaussian_prior_name):
        scale = data.mcmc_parameters[nuisance_name]['scale']
        chi2 += (data.mcmc_parameters[nuisance_name]['current'] * scale - self.gaussian_prior_center[idx_nuisance])**2 / self.gaussian_prior_sigma[idx_nuisance]**2

return -chi2 / 2.

Cheers, Fabian