Made a new branch called logLmax with the updated anesthetic.examples…

[DilyOng]

….perfect_ns.correlated_gaussian.

Description

A bug was found in the function anesthetic.examples.perfect_ns.correlated_gaussian(). The generated Gaussian likelihood is not normalised and the evidence is not unity. Three changes have been made.

• The function description is updated to remind the users that the correlated_gaussian() does not produce a correlated gaussian likelihood which is normalised and the evidence is not unity.
• An extra argument called LogLmax is introduced to the function correlated_gaussian().
• Within the correlated_gaussian() function, the definition of logLike() is corrected.

Fixes # (issue)

Checklist:

• I have performed a self-review of my own code
• My code is PEP8 compliant (flake8 anesthetic tests)
• My code contains compliant docstrings (pydocstyle --convention=numpy anesthetic)
• New and existing unit tests pass locally with my changes (python -m pytest)
• I have added tests that prove my fix is effective or that my feature works
• I have appropriately incremented the semantic version number in both README.rst and anesthetic/_version.py

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (7c30e72) 100.00% compared to head (44f6dd5) 100.00%.

Additional details and impacted files

@@            Coverage Diff            @@
##            master      #348   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           34        34           
  Lines         2984      2978    -6     
=========================================
- Hits          2984      2978    -6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[williamjameshandley]

Many thanks for this proposal @DilyOng.

At the moment the tests are failing because your changes are not flake8 PEP compliant. You can see this by clicking 'Details' next to the failed test, and see it locally by running python -m flake8 anesthetic tests. You also need to increment the version number as the Checklist at the top of the PR suggests.

Could you add a test to tests.test_examples.test_perfect_ns_correlated_gaussian which checks that the logLmax functionality is working. The best thing to do would be to make a run with logLmax equal to some number (e.g. 10), and check that samples.logL.mean() is close to logLmax-d/2, where d is the number of parameters.

One other thing we should consider is putting in a logLmax into the gaussian function as well. This would also need a test.

[williamjameshandley]

Many thanks @DilyOng.

1. Could you remove the tests/tests.ipynb file
2. Could you move test_perfect_ns_correlated_gaussian.py int othe file tests/test_examples.py
3. Could you add the logLmax variable to gaussian
4. The version number should be 2.5.0 as we have added functionality
5. Could you modify & simplify the planck_gaussian function so that it uses your new functionality!

[williamjameshandley]

Hi @DilyOng,

I have pushed a minor but structural change, which implements the functionality that in gaussian, wedding_cake and correlated_gaussian the 'rest of the arguments' *args, **kwargs are passed on to the NestedSamples constructor.

The main reason for doing this is that it allows us to pass columns and params to correlated_gaussian in the planck_gaussian function.

One further minor change is that with logL_max being set means that the logL_mean comes with an additional nested sampling error, so the tests associated with logL_mean` had to be adjusted.

If you're happy with these changes, you should add yourself to the author list in pyproject.toml, and then I'll approve these changes to be merged.

[williamjameshandley]

Many thanks @DilyOng. Please press 'squash and merge'.

[lukashergt]

Hi @DilyOng, just one request from me. Could we add a docstring entries for logLmax, please?

[lukashergt]

logLmax is lacking a docstring entry.

[DilyOng]

Thanks for your comment. But logLmax needs to be a number. I'm not sure what you mean I am sorry...

[lukashergt]

The docstring underneath the function definition is the text that describes what the function does such that a new user can quickly inform themselves on the correct usage. Such a function docstring should include a list of (at least) all the explicit function arguments. By numpy convention (which is the one that we are using), all these function arguments get listed underneath a "Parameters" heading. In jupyter notebooks this docstring is easily accessible by pressing shift+tab inside a function.

Here, the docstring currently is:

    """Perfect nested sampling run for a spherical Gaussian & prior.

    Up to normalisation this is identical to the example in John Skilling's
    original paper https://doi.org/10.1214/06-BA127 . Both spherical Gaussian
    width sigma and spherical uniform prior width R are centered on zero

    Parameters
    ----------
    nlive : int
        number of live points#

    ndims : int
        dimensionality of gaussian

    sigma : float
        width of gaussian likelihood

    R : float
        radius of gaussian prior

    logLmin : float
        loglikelihood at which to terminate

    The remaining arguments are passed to the
    :class:`anesthetic.samples.NestedSamples` constructor.

    Returns
    -------
    samples : :class:`anesthetic.samples.NestedSamples`
        Nested sampling run
    """

So this currently lists nlive, ndims, sigma, R, and logLmin, but it lacks a description of logLmax.

[DilyOng]

Thank you very much for explaining! I've implemented your suggestion.

[lukashergt]

Hmm, I was hoping for a bit more description...

Then again, the docstring for the other arguments is also very minimalistic... and possibly a bit misleading... Is this really a nested sampling run of a spherical Gaussian and prior? Or would it be more accurate to call it a nested sampling run of a spherical Gaussian likelihood from a spherical top-hat prior? In which case the docstring for R should be rather the radius of the spherical top-hat prior...

Anyhow, it's not right to ask you to clean all that up (unless you want to), but it would be nice if you could make the logLmax docstring a bit more descriptive. Something along the lines that, yes, this defines the maximum of the Gaussian likelihood, but in effect this only shifts the log-likelihood up and down by a constant. Is there a reason anyone should choose a value different from the default logLmax=0? What prompted you to introduce this argument the first place? If that is a gotcha that would be good for other people to be aware of, then this is a good place to warn about it.

[DilyOng]

Hi @lukashergt I am terribly sorry for my late reply, I somehow missed your newest comment.

The reason why we need to add LogLmax is because we’ve found an issue with the definition of combined Gaussian likelihood. We are recently adding a new function (in the branch called tension) to anesthetic for calculating tension statistics between datasets. The initial test involves creating 2 mock datasets A and B with d-dimensional Gaussian likelihood and a combined dataset AB (combining A and B on the likelihood level). We discovered that the quantity of logLmax is actually needed explicitly in the calculation of combined likelihoods. Besides, we have noticed that the correlated_gaussian() function does not produce a normalised correlated gaussian likelihood and the evidence isn’t unity. Therefore, we need to update the function description. Please refer to equations 14-21 in the paper Quantifying tensions in cosmological parameters: Interpreting the DES evidence ratio [1902.04029]. Thank you very much!

[lukashergt]

docstring entry for logLmax

[DilyOng]

Done!

Please implement changes requested by @!ukashergt

[lukashergt]

Sorry, I did not mean to halt this completely, just thought it best to document things thoroughly while the reason for changes/additions are fresh in our minds, but not at the cost of causing a standstill, so I'd rather approve this now.