Prior for derived parameters [of theory code] not working

[vivianmiranda]

Hi Again

I was trying to reproduce Fig 7 of

https://arxiv.org/pdf/1806.04649.pdf

So I created a fake likelihood that just return chi^2 = 0.0 and I made the following yaml

test6.txt

which according to instruction (from H0) should have limited the range of omega_m

omegam: max: 1 min: 0 latex: '\Omega_\mathrm{m}'

but when I see the chain files - I see that there are bunch of accepted points with omegam outside this range.

Here one example

TEST_PRIORS_1.txt

I can do similar thing in the likelihood (return -\infty if omegam is outside the range) but it seems awkward that min/max does not work (or only works for H0)

Best

[JesusTorrado]

Hi Vivian,

Derived parameters do not have priors as such: max and min are just limits passed to GetDist for plotting and computing c.l. intervals.

If you want to impose this as a prior, you have to define the equivalent prior in the prior block. Normally, you could use omega_m in the prior directly and make something like lambda omega_m: 0 < omega_m < 1.

Unfortunately, in this particular case of yours, since omega_mis not an input/sampled parameter, but is instead requested to the theory code, and that one has not yet run at prior evaluation, you need to derived it from input parameters, so it would be lambda ombh2, omch2, mmu: [don't remember now].

Working around that limitation is on my list. I'll try to get it done today. EDIT: a little harder than I expected. Need to think more about it

[JesusTorrado]

Another alternative: add the following likelihood (likelihoods are evaluated after the theory code)

omb: "lambda _theory={'omegam': None}: np.log(0 <= _theory.get_param('omegam') <= 1)"

[vivianmiranda]

Thank you!

It would be nice to have these priors before the theory code - like in CosmoMC.

For example - there are cases in model extensions where CAMB could break down or be extremely slow if some derived parameters are this or that.

In my inflation paper with Wayne - the calculation was just wrong if something happened with a combination of parameters. Having them to be evaluated before theory saves loads of time.

Thank you again for all support.

PS: So why the H0 min max work?

Best Vivian

On Dec 17, 2018, at 4:52 AM, Jesús Torrado <[email protected]mailto:[email protected]> wrote:

Another alternative: add the following likelihood (likelihoods are evaluated after the theory code)

omb: "lambda _theory={'omegam': None}: np.log(0 <= _theory.get_param('omegam') <= 1)"

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/CobayaSampler/cobaya/issues/18#issuecomment-447818811, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADD96KKBSQ14nqgKpIJm5CTWv4qgVu1jks5u54V5gaJpZM4ZUdUm.

[JesusTorrado]

Hi Vivian,

In order to have them evaluated before the theory code, Cobaya needs to know how to extract them from the input and sampled parameters. To do that, you need to define a prior as in my first suggestion.

An alternative would be for Cobaya to implement those definitions internally, as CosmoMC does. I am not keen at all on that approach (which is what CosmoMC does), since I strongly prefer to keep Cobaya as a "general sampler with Cosmology wrappers", agnostic to parameter names. Otherwise, this opens a can of worms: what's your Omega_m definition? Are there extra neutrino species or any other contribution? If I have to write the full logic into Cobaya, I may end up rewriting the "background" part of CAMB/CLASS, which is way harder to maintain, and not the job of a sampler anyway. It would make more sense to turn CAMB/CLASS into a library whose parts are called sequentially and only on demand, e.g. background quantities only as long as no perturbations-related quantity is requested. (I hope that did not sound as I am annoyed by the suggestion; just making clear the logic behind my design decisions, with which Antony does not always agree BTW).

In any case, if you get Omega_m > 1 without any of the workarounds above but you imposed 0-curvature (the default), this is actually a CAMB issue.

The H0 case appears to work because there are equivalent limitations inside CAMB, which are the actual ones being enforced. If you made the min/max larger, it should have no effect.

[JussiValiviita]

If you want to impose this as a prior, you have to define the equivalent prior in the prior block. Normally, you could use omega_m in the prior directly and make something like lambda omega_m: 0 < omega_m < 1.

Dear Jesus, I have a bit simpler problem than Vivian. My derived parameter (e.g., 'we' for the early-time dark energy equation of state) is calculated directly from the primary parameters (e.g., w and wa) without a need for a camb call: we = w + wa. Then I would like to implement a "top hat" prior on this derived parameter; for example, -3 < we < -0.34. However, I have not found any way to do this (without touching the actual python and fortran codes), though I think this kind of prior on a linear combination of two primary parameters is a so commonly needed trick that it should be doable directly by using the yaml files. (I did not find any examples in the documentation, but please accept my apology and point me to the right subpage if an example exists!)

Among others, I have tried the following in the main yaml file of my run: prior: 'lambda w,wa: -3 < w+wa < -0.34' In most places of yaml this line causes a crash, but if I put it in the "likelihood" block (of the main yaml) with a right indentation, the run starts, but very soon I start to see in the created chain files values that are we > -0.34, i.e., the prior is ignored.

So, what do you mean by a "prior block" in the quoted comment? Where exactly should I place my prior?

And to clarify, I have a couple of similar priors in my model, and they need to be implemented BEFORE the camb call since they are mostly used to prevent unphysical combinations of parameters that will cause a segmentation fault or other unrecoverable crash of camb. Thus implementing the prior in the likelihood that is called AFTER camb is not possible. The prior should be checked immediately once the primary parameters have been drawn. How to do this?

Thanks in advance for your time, Jussi Väliviita

[cmbant]

You can put a prior block at the highest level, e.g. see how the Planck likleihoods use prior_SZ_2015. The block is a list of named priors.

[fbianchini]

Hi there, I'm encountering a similar problem. Specifically, I'm trying to add an external prior on omegam but Jesus's suggestion to include an additional "likelihood" as

like_omm: "lambda _theory={'omegam': None}: stats.norm.logpdf(_theory.get_param('omegam'), loc=0.334, scale=0.018)"

doesn't seem to do the trick as it spits out the following error:

[model] *ERROR* Component 'like_omm' seems not to depend on any parameters (neither directly nor indirectly)

For completeness, I'm using CAMB as theory code and defining omegam in the params block:

omegam:
latex: \Omega_\mathrm{m}

Did something change in more recent versions or am I incorrectly adding the prior?

Thanks for any and all inputs!

Cheers, ~federico

[cmbant]

I think the syntax is now (from memory) something like

like_omm:  
  external: "lambda _self: stats.norm.logpdf(_self.get_param('omegam'), loc=0.334, scale=0.018)"
  requires: omegam

[fbianchini]

Thanks for the prompt reply, Antony!

I've updated the syntax accordingly but the code is now complaining that it cannot retrieve the calculated quantities (as shown in the error message below)

[like_omm] Got parameters {'omegam': 0.3347663875519032}
[like_omm] Computing new state
[like_omm] *ERROR* Cannot retrieve calculated quantities: nothing has been computed yet (maybe the prior was -infinity?)
[like_omm] External function failed at evaluation.

Any pointers would be greatly appreciated!

~federico

[cmbant]

Not sure, can you attach a complete example to reproduce?

[fbianchini]

For sure, here's a link to the yaml file I've been using, I'm sampling over a 5 param LCDM cosmology (no tau) including the Planck PR4 lensing dataset + a prior on Omega_M.

[cmbant]

Sorry, I think syntax should be _self.provider.get_param()

[fbianchini]

Thanks heaps Antony, that worked like a charm!