Tension estimate for non-Gaussian distributions¶

Cyrille Doux (doux@lpsc.in2p3.fr), Marco Raveri (marco.raveri@unige.it)

In this notebook we show how to calculate the level of tension between two experiments.

In particular we show how to compute the statistical significance of a parameter shift between two experiments, DES Y1 and Planck 18, with the two techniques discussed in Raveri and Doux (2021), arXiv:2105.03324.

In [1]:
# Show plots inline, and load main getdist plot module and samples class
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext autoreload
%autoreload 1
# import libraries:
import sys, os
os.environ['TF_USE_LEGACY_KERAS'] = '1'  # needed for tensorflow KERAS compatibility
os.environ['DISPLAY'] = 'inline'  # hack to get getdist working
sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'../..')))
from getdist import plots
import getdist
getdist.chains.print_load_details = False
import matplotlib.pyplot as plt

# import the tensiometer tools that we need:
import tensiometer
from tensiometer.utilities import stats_utilities as utilities
from tensiometer import mcmc_tension

We start importing the relevant samples. The code requires the two chains to be getdist MCSamples objects. Check GetDist documentation for example codes of how to import or create them.

In [2]:
# load the samples (remove no burn in since the example chains have already been cleaned):
chains_dir = os.path.realpath(os.path.join(os.getcwd(), '../..', 'test_chains'))
# the Planck 2018 TTTEEE chain:
chain_1 = getdist.mcsamples.loadMCSamples(file_root=os.path.join(chains_dir, 'Planck18TTTEEE'), no_cache=True)
# the DES Y1 3x2 chain:
chain_2 = getdist.mcsamples.loadMCSamples(file_root=os.path.join(chains_dir, 'DES'), no_cache=True)

Now we compute the parameter difference distribution.

If we denote $P_1(\theta_1)$ and $P_2(\theta_2)$ as the two independent distributions, then, defining $\Delta\theta \equiv \theta_1 - \theta_2$ then the distribution of parameter differences, $\Delta\theta$ is given by:

\begin{align} P(\Delta \theta) = \int P_1(\theta)P_2(\theta - \Delta\theta) \, d\theta \end{align}

Samples from this distribution can be computed as differences between samples of the two distributions. This usually results in a large number of samples so, by default, we undersample. The boost parameter increases the samples retained, up to $O(n^2)$, and should be set to ensure that we have enough samples for later.

In [3]:
# the difference chain:
diff_chain = mcmc_tension.parameter_diff_chain(chain_1, chain_2, boost=4)

Now we do a sanity check plot:

In [4]:
param_names = diff_chain.getParamNames().getRunningNames()
g = plots.get_subplot_plotter()
g.triangle_plot([diff_chain], params=param_names, filled=True, markers={_p:0 for _p in param_names})
plt.show()
No description has been provided for this image

Looks non-Gaussian uh? Let's see below how to cope with this!

We need to calculate this integral:

\begin{align} \Delta \equiv \int_{P(\Delta\theta) > P(0)} P(\Delta\theta) \, d\Delta \theta \end{align}

which gives the probability mass enclosed in the (full-D) iso-contour that touches the value corresponding to zero shift.

The main problem of this integral is that we have samples from the parameter difference distribution but we cannot associate a probability value to those samples. There are two ways to cope with this problem, that we outline below.

Normalizing flow estimate of parameter shifts:¶

The first strategy is to build a normalizing flow model for the parameter difference distribution. We learn from the samples the mapping between parameter spaces that Gaussianizes the distribution. Once this is done we can sample from the learned distribution and compute probability values. The previous integral is then Monte-Carlo integrated.

The code provides a helper function, tensiometer.mcmc_tension.flow_parameter_shift(diff_chain), to create the model, train it and compute the shift significance. We show here how to use it and later how to proceed manually.

Training the normalizing flow is going to take a little while. This is a good time to take break!

In [5]:
# the helper function passes all keyword arguments downstream so we can easily set manually all options
kwargs = {
          'feedback': 2,  # verbosity level
          'pop_size': 1,  # controls the number of flows that are trained. The best one is then selected. Since time to solution scales linearly with pop_size, it is a good idea to set it to 1 for a first run. 
        }

# call the helper function:
results, diff_flow = tensiometer.mcmc_tension.flow_parameter_shift(diff_chain,   # parameter difference chain
                                      cache_dir=None,  # directory where the trained flow is saved, to save time next time
                                      root_name='sprob',  # name of the cached flow (usefull to have several in the same folder)
                                      **kwargs
                                      )

# unpack results:
shift_P, shift_hi, shift_low = results
No description has been provided for this image 
20/20 - 9s - loss: 6.9615 - val_loss: 6.9724 - lr: 0.0010 - 9s/epoch - 463ms/step
* Population optimizer:
    - best model is number 1
    - best loss function is 6.96
    - best validation loss function is 6.97
    - population losses [6.97]
In [6]:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99153 +-0.00030 --0.00029
    n_sigma = 2.633 +-0.012 --0.012

Manual flow calculation¶

Now we see how to unpack the previous helper function and calculate a tension from scratch.

We first create and train the synthetic probability model. Again this is going to take a bit. Normally (not in a tutorial) one would save the flow model to cache to avoid retraining every time...

In [7]:
kwargs = {
          'feedback': 1,
          'plot_every': 1000,
          'pop_size': 1,
        }

flow = tensiometer.synthetic_probability.synthetic_probability.flow_from_chain(diff_chain,  # parameter difference chain
                                                         **kwargs)

* Initializing samples
    - time taken: 0.0382 seconds
* Initializing fixed bijector
    - time taken: 0.3564 seconds
* Initializing trainable bijector
    - time taken: 1.5621 seconds
* Initializing training dataset
    - time taken: 2.0113 seconds
* Initializing transformed distribution
    - time taken: 0.0130 seconds
* Initializing loss function
    - time taken: 0.0000 seconds
* Initializing training model
    - Compiling model
    - time taken: 0.0592 seconds
    - time taken: 1.2308 seconds
* Training
    - Compiling model
    - time taken: 0.0527 seconds
Epoch 1/100
20/20 - 22s - loss: 8.4430 - val_loss: 8.3568 - lr: 0.0010 - 22s/epoch - 1s/step
Epoch 2/100
20/20 - 6s - loss: 8.2944 - val_loss: 8.2207 - lr: 0.0010 - 6s/epoch - 294ms/step
Epoch 3/100
20/20 - 6s - loss: 8.1576 - val_loss: 8.0842 - lr: 0.0010 - 6s/epoch - 310ms/step
Epoch 4/100
20/20 - 6s - loss: 8.0343 - val_loss: 7.9821 - lr: 0.0010 - 6s/epoch - 307ms/step
Epoch 5/100
20/20 - 6s - loss: 7.9486 - val_loss: 7.9045 - lr: 0.0010 - 6s/epoch - 294ms/step
Epoch 6/100
20/20 - 6s - loss: 7.8745 - val_loss: 7.8401 - lr: 0.0010 - 6s/epoch - 289ms/step
Epoch 7/100
20/20 - 6s - loss: 7.8234 - val_loss: 7.7970 - lr: 0.0010 - 6s/epoch - 292ms/step
Epoch 8/100
20/20 - 6s - loss: 7.7855 - val_loss: 7.7676 - lr: 0.0010 - 6s/epoch - 319ms/step
Epoch 9/100
20/20 - 7s - loss: 7.7528 - val_loss: 7.7342 - lr: 0.0010 - 7s/epoch - 340ms/step
Epoch 10/100
20/20 - 6s - loss: 7.7211 - val_loss: 7.7027 - lr: 0.0010 - 6s/epoch - 300ms/step
Epoch 11/100
20/20 - 6s - loss: 7.6888 - val_loss: 7.6714 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 12/100
20/20 - 6s - loss: 7.6587 - val_loss: 7.6429 - lr: 0.0010 - 6s/epoch - 303ms/step
Epoch 13/100
20/20 - 7s - loss: 7.6323 - val_loss: 7.6182 - lr: 0.0010 - 7s/epoch - 330ms/step
Epoch 14/100
20/20 - 6s - loss: 7.6089 - val_loss: 7.5970 - lr: 0.0010 - 6s/epoch - 311ms/step
Epoch 15/100
20/20 - 6s - loss: 7.5884 - val_loss: 7.5769 - lr: 0.0010 - 6s/epoch - 293ms/step
Epoch 16/100
20/20 - 6s - loss: 7.5681 - val_loss: 7.5562 - lr: 0.0010 - 6s/epoch - 321ms/step
Epoch 17/100
20/20 - 6s - loss: 7.5460 - val_loss: 7.5324 - lr: 0.0010 - 6s/epoch - 299ms/step
Epoch 18/100
20/20 - 6s - loss: 7.5191 - val_loss: 7.5013 - lr: 0.0010 - 6s/epoch - 307ms/step
Epoch 19/100
20/20 - 6s - loss: 7.4826 - val_loss: 7.4585 - lr: 0.0010 - 6s/epoch - 318ms/step
Epoch 20/100
20/20 - 6s - loss: 7.4341 - val_loss: 7.4029 - lr: 0.0010 - 6s/epoch - 308ms/step
Epoch 21/100
20/20 - 6s - loss: 7.3758 - val_loss: 7.3441 - lr: 0.0010 - 6s/epoch - 324ms/step
Epoch 22/100
20/20 - 6s - loss: 7.3205 - val_loss: 7.2937 - lr: 0.0010 - 6s/epoch - 283ms/step
Epoch 23/100
20/20 - 6s - loss: 7.2736 - val_loss: 7.2512 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 24/100
20/20 - 6s - loss: 7.2341 - val_loss: 7.2144 - lr: 0.0010 - 6s/epoch - 302ms/step
Epoch 25/100
20/20 - 7s - loss: 7.1988 - val_loss: 7.1805 - lr: 0.0010 - 7s/epoch - 326ms/step
Epoch 26/100
20/20 - 7s - loss: 7.1665 - val_loss: 7.1492 - lr: 0.0010 - 7s/epoch - 326ms/step
Epoch 27/100
20/20 - 6s - loss: 7.1370 - val_loss: 7.1203 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 28/100
20/20 - 6s - loss: 7.1105 - val_loss: 7.0977 - lr: 0.0010 - 6s/epoch - 288ms/step
Epoch 29/100
20/20 - 6s - loss: 7.0892 - val_loss: 7.0761 - lr: 0.0010 - 6s/epoch - 310ms/step
Epoch 30/100
20/20 - 6s - loss: 7.0728 - val_loss: 7.0614 - lr: 0.0010 - 6s/epoch - 307ms/step
Epoch 31/100
20/20 - 6s - loss: 7.0595 - val_loss: 7.0504 - lr: 0.0010 - 6s/epoch - 323ms/step
Epoch 32/100
20/20 - 6s - loss: 7.0501 - val_loss: 7.0428 - lr: 0.0010 - 6s/epoch - 314ms/step
Epoch 33/100
20/20 - 6s - loss: 7.0420 - val_loss: 7.0349 - lr: 0.0010 - 6s/epoch - 307ms/step
Epoch 34/100
20/20 - 6s - loss: 7.0366 - val_loss: 7.0302 - lr: 0.0010 - 6s/epoch - 293ms/step
Epoch 35/100
20/20 - 6s - loss: 7.0309 - val_loss: 7.0272 - lr: 0.0010 - 6s/epoch - 306ms/step
Epoch 36/100
20/20 - 6s - loss: 7.0266 - val_loss: 7.0204 - lr: 0.0010 - 6s/epoch - 310ms/step
Epoch 37/100
20/20 - 6s - loss: 7.0217 - val_loss: 7.0167 - lr: 0.0010 - 6s/epoch - 290ms/step
Epoch 38/100
20/20 - 6s - loss: 7.0179 - val_loss: 7.0126 - lr: 0.0010 - 6s/epoch - 299ms/step
Epoch 39/100
20/20 - 6s - loss: 7.0144 - val_loss: 7.0097 - lr: 0.0010 - 6s/epoch - 311ms/step
Epoch 40/100
20/20 - 6s - loss: 7.0107 - val_loss: 7.0062 - lr: 0.0010 - 6s/epoch - 313ms/step
Epoch 41/100
20/20 - 6s - loss: 7.0077 - val_loss: 7.0034 - lr: 0.0010 - 6s/epoch - 304ms/step
Epoch 42/100
20/20 - 6s - loss: 7.0050 - val_loss: 7.0013 - lr: 0.0010 - 6s/epoch - 314ms/step
Epoch 43/100
20/20 - 8s - loss: 7.0025 - val_loss: 6.9980 - lr: 0.0010 - 8s/epoch - 375ms/step
Epoch 44/100
20/20 - 6s - loss: 6.9996 - val_loss: 6.9957 - lr: 0.0010 - 6s/epoch - 306ms/step
Epoch 45/100
20/20 - 6s - loss: 6.9977 - val_loss: 6.9931 - lr: 0.0010 - 6s/epoch - 304ms/step
Epoch 46/100
20/20 - 8s - loss: 6.9960 - val_loss: 6.9910 - lr: 0.0010 - 8s/epoch - 389ms/step
Epoch 47/100
20/20 - 8s - loss: 6.9937 - val_loss: 6.9892 - lr: 0.0010 - 8s/epoch - 406ms/step
Epoch 48/100
20/20 - 7s - loss: 6.9913 - val_loss: 6.9878 - lr: 0.0010 - 7s/epoch - 339ms/step
Epoch 49/100
20/20 - 6s - loss: 6.9896 - val_loss: 6.9854 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 50/100
20/20 - 6s - loss: 6.9883 - val_loss: 6.9858 - lr: 0.0010 - 6s/epoch - 324ms/step
Epoch 51/100
20/20 - 6s - loss: 6.9872 - val_loss: 6.9827 - lr: 0.0010 - 6s/epoch - 283ms/step
Epoch 52/100
20/20 - 6s - loss: 6.9848 - val_loss: 6.9809 - lr: 0.0010 - 6s/epoch - 289ms/step
Epoch 53/100
20/20 - 6s - loss: 6.9833 - val_loss: 6.9793 - lr: 0.0010 - 6s/epoch - 310ms/step
Epoch 54/100
20/20 - 6s - loss: 6.9820 - val_loss: 6.9802 - lr: 0.0010 - 6s/epoch - 302ms/step
Epoch 55/100
20/20 - 6s - loss: 6.9807 - val_loss: 6.9768 - lr: 0.0010 - 6s/epoch - 299ms/step
Epoch 56/100
20/20 - 6s - loss: 6.9796 - val_loss: 6.9767 - lr: 0.0010 - 6s/epoch - 319ms/step
Epoch 57/100
20/20 - 6s - loss: 6.9786 - val_loss: 6.9740 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 58/100
20/20 - 6s - loss: 6.9767 - val_loss: 6.9732 - lr: 0.0010 - 6s/epoch - 317ms/step
Epoch 59/100
20/20 - 6s - loss: 6.9760 - val_loss: 6.9728 - lr: 0.0010 - 6s/epoch - 311ms/step
Epoch 60/100
20/20 - 6s - loss: 6.9747 - val_loss: 6.9715 - lr: 0.0010 - 6s/epoch - 311ms/step
Epoch 61/100
20/20 - 6s - loss: 6.9733 - val_loss: 6.9704 - lr: 0.0010 - 6s/epoch - 309ms/step
Epoch 62/100
20/20 - 6s - loss: 6.9722 - val_loss: 6.9685 - lr: 0.0010 - 6s/epoch - 316ms/step
Epoch 63/100
20/20 - 6s - loss: 6.9711 - val_loss: 6.9684 - lr: 0.0010 - 6s/epoch - 320ms/step
Epoch 64/100
20/20 - 6s - loss: 6.9704 - val_loss: 6.9663 - lr: 0.0010 - 6s/epoch - 290ms/step
Epoch 65/100
20/20 - 6s - loss: 6.9691 - val_loss: 6.9656 - lr: 0.0010 - 6s/epoch - 290ms/step
Epoch 66/100
20/20 - 6s - loss: 6.9676 - val_loss: 6.9641 - lr: 0.0010 - 6s/epoch - 299ms/step
Epoch 67/100
20/20 - 7s - loss: 6.9668 - val_loss: 6.9645 - lr: 0.0010 - 7s/epoch - 328ms/step
Epoch 68/100
20/20 - 6s - loss: 6.9662 - val_loss: 6.9624 - lr: 0.0010 - 6s/epoch - 309ms/step
Epoch 69/100
20/20 - 6s - loss: 6.9648 - val_loss: 6.9627 - lr: 0.0010 - 6s/epoch - 305ms/step
Epoch 70/100
20/20 - 6s - loss: 6.9645 - val_loss: 6.9601 - lr: 0.0010 - 6s/epoch - 302ms/step
Epoch 71/100
20/20 - 6s - loss: 6.9637 - val_loss: 6.9591 - lr: 0.0010 - 6s/epoch - 297ms/step
Epoch 72/100
20/20 - 6s - loss: 6.9618 - val_loss: 6.9584 - lr: 0.0010 - 6s/epoch - 311ms/step
Epoch 73/100
20/20 - 6s - loss: 6.9610 - val_loss: 6.9580 - lr: 0.0010 - 6s/epoch - 296ms/step
Epoch 74/100
20/20 - 6s - loss: 6.9594 - val_loss: 6.9561 - lr: 0.0010 - 6s/epoch - 294ms/step
Epoch 75/100
20/20 - 6s - loss: 6.9590 - val_loss: 6.9570 - lr: 0.0010 - 6s/epoch - 314ms/step
Epoch 76/100
20/20 - 7s - loss: 6.9583 - val_loss: 6.9575 - lr: 0.0010 - 7s/epoch - 328ms/step
Epoch 77/100
20/20 - 7s - loss: 6.9574 - val_loss: 6.9538 - lr: 0.0010 - 7s/epoch - 366ms/step
Epoch 78/100
20/20 - 8s - loss: 6.9557 - val_loss: 6.9547 - lr: 0.0010 - 8s/epoch - 392ms/step
Epoch 79/100
20/20 - 8s - loss: 6.9562 - val_loss: 6.9526 - lr: 0.0010 - 8s/epoch - 414ms/step
Epoch 80/100
20/20 - 7s - loss: 6.9541 - val_loss: 6.9508 - lr: 0.0010 - 7s/epoch - 371ms/step
Epoch 81/100
20/20 - 7s - loss: 6.9539 - val_loss: 6.9499 - lr: 0.0010 - 7s/epoch - 358ms/step
Epoch 82/100
20/20 - 7s - loss: 6.9521 - val_loss: 6.9494 - lr: 0.0010 - 7s/epoch - 368ms/step
Epoch 83/100
20/20 - 7s - loss: 6.9512 - val_loss: 6.9479 - lr: 0.0010 - 7s/epoch - 355ms/step
Epoch 84/100
20/20 - 7s - loss: 6.9499 - val_loss: 6.9473 - lr: 0.0010 - 7s/epoch - 349ms/step
Epoch 85/100
20/20 - 7s - loss: 6.9497 - val_loss: 6.9473 - lr: 0.0010 - 7s/epoch - 348ms/step
Epoch 86/100
20/20 - 7s - loss: 6.9490 - val_loss: 6.9476 - lr: 0.0010 - 7s/epoch - 361ms/step
Epoch 87/100
20/20 - 8s - loss: 6.9483 - val_loss: 6.9460 - lr: 0.0010 - 8s/epoch - 392ms/step
Epoch 88/100
20/20 - 7s - loss: 6.9480 - val_loss: 6.9442 - lr: 0.0010 - 7s/epoch - 356ms/step
Epoch 89/100
20/20 - 8s - loss: 6.9466 - val_loss: 6.9438 - lr: 0.0010 - 8s/epoch - 377ms/step
Epoch 90/100
20/20 - 7s - loss: 6.9449 - val_loss: 6.9430 - lr: 0.0010 - 7s/epoch - 353ms/step
Epoch 91/100
20/20 - 7s - loss: 6.9451 - val_loss: 6.9422 - lr: 0.0010 - 7s/epoch - 365ms/step
Epoch 92/100
20/20 - 8s - loss: 6.9438 - val_loss: 6.9407 - lr: 0.0010 - 8s/epoch - 388ms/step
Epoch 93/100
20/20 - 8s - loss: 6.9433 - val_loss: 6.9403 - lr: 0.0010 - 8s/epoch - 402ms/step
Epoch 94/100
20/20 - 7s - loss: 6.9422 - val_loss: 6.9393 - lr: 0.0010 - 7s/epoch - 371ms/step
Epoch 95/100
20/20 - 7s - loss: 6.9416 - val_loss: 6.9387 - lr: 0.0010 - 7s/epoch - 367ms/step
Epoch 96/100
20/20 - 7s - loss: 6.9405 - val_loss: 6.9382 - lr: 0.0010 - 7s/epoch - 375ms/step
Epoch 97/100
20/20 - 7s - loss: 6.9398 - val_loss: 6.9372 - lr: 0.0010 - 7s/epoch - 354ms/step
Epoch 98/100
20/20 - 8s - loss: 6.9392 - val_loss: 6.9363 - lr: 0.0010 - 8s/epoch - 404ms/step
Epoch 99/100
20/20 - 7s - loss: 6.9394 - val_loss: 6.9374 - lr: 0.0010 - 7s/epoch - 367ms/step
Epoch 100/100
20/20 - 6s - loss: 6.9383 - val_loss: 6.9359 - lr: 0.0010 - 6s/epoch - 325ms/step
In [8]:
# we can plot training summaries to make sure training went smoothly:
flow.training_plot()
plt.show();
No description has been provided for this image
In [9]:
# calculate the shift probability:
shift_P, shift_low, shift_hi = mcmc_tension.estimate_shift(flow)
# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99655 +0.00018 -0.00019
    n_sigma = 2.925 +0.017 -0.017

Note that there is some variance in the result due to initialization. You can repeat the above calculation some times to evaluate the variance.

We can now plot the learned distribution. To do so we draw some samples from the learned distribution and then feed them to getdist plotting:

In [10]:
# sample the flow:
N = 10000
flow_samples = flow.MCSamples(N)

# build Gaussian approximation:
gaussian_approx = tensiometer.gaussian_tension.gaussian_approximation(diff_chain)
In [11]:
colors=['orange', 'dodgerblue', 'k']
g = plots.get_subplot_plotter()
g.settings.num_plot_contours = 2
g.triangle_plot([diff_chain, flow_samples, gaussian_approx], params=param_names,
                filled=False, markers={_p:0 for _p in param_names},
                colors=colors, diag1d_kwargs={'colors':colors})
plt.show();
No description has been provided for this image

As we can see the two distributions match astonishingly well, capturing all non-Gaussian feature that we can identify in this plot.

If you are interested in more details about normalizing flow modeling of probability distributions check out the synthetic probability tutorial!

KDE estimate of parameter shifts:¶

We first run the KDE algorithm with default settings and high feedback to have a sense of its inner workings:

In [12]:
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift(diff_chain, feedback=10)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Dimension       : 6
N    samples    : 1339087
Neff samples    : 1147971.95
Smoothing scale : [0.056 0.056 0.056 0.056 0.056 0.056]
Building KD-Tree with leafsize = 400
Neighbours elimination
neighbor_elimination: chunk 1
    surviving elements 5540 of 1257016
neighbor_elimination: chunk 2
    surviving elements 4896 of 1257016
neighbor_elimination: chunk 3
    surviving elements 4699 of 1257016
neighbor_elimination: chunk 4
    surviving elements 4626 of 1257016
neighbor_elimination: chunk 5
    surviving elements 4589 of 1257016
neighbor_elimination: polishing

100%|##########| 4589/4589 [02:26<00:00, 31.30it/s]

    surviving elements 4495 of 1257016
KDE method: neighbor_elimination
Time taken for KDE calculation: 238.2 (s)
Shift probability considering all parameters:
    P = 0.99659 +0.00005 -0.00005
    n_sigma = 2.928 +0.005 -0.005

As we can see the algorithm proceeds in two steps, elimination of points that are clearly above the probability of zero shift based on the probability estimate containing just a few nearest points and then brute force polishing of the leftover points. This gets the right answer and takes a fairly reasonable amount of time.

You could try to run this cell with method='brute_force' to see the type of performance improvement that this algorithm achieves.

One key parameter in the estimate of the KDE shift is the smoothing scale for the pdf. The default choice is the one that minimizes the mean integrated square error (MISE) under assumptions of Gaussianity.

In [13]:
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift(diff_chain, scale='MISE', feedback=0)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99659 +0.00005 -0.00005
    n_sigma = 2.928 +0.005 -0.005

But we can try different choices. In particular we could try the asyntotic MISE estimator (AMISE):

In [14]:
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift(diff_chain, scale='AMISE', feedback=0)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99738 +0.00004 -0.00004
    n_sigma = 3.009 +0.005 -0.005

which is slightly undersmoothing and hence resulting in a slightly higher tension.

Or we could try the maximum bandwidth that is (by design) oversmoothing:

In [15]:
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift(diff_chain, scale='MAX', feedback=0)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99181 +0.00008 -0.00008
    n_sigma = 2.644 +0.003 -0.003

As we can see the MISE smoothing scale achieves a balance between these two. Overestimating the smoothing scale usually results in slightly smaller tensions.

We can now try the adaptive bandwidth:

In [16]:
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift(diff_chain, scale='BALL', feedback=0)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

Shift probability considering all parameters:
    P = 0.99593 +0.00005 -0.00006
    n_sigma = 2.873 +0.004 -0.004

As we can see all results are mostly in agreement with each other and this provides an invaluable cross check of the different calculation techniques.

The number of samples used in the calculation can be changed to make sure that everything is converged and makes sense. Since all algorithms work in a fairly straightforward way on a laptop we encourage to always try and compare different outputs to make sure the result is sensible.

We now verify that the calculation makes physical sense. We know that a difference between these two results mostly lives in a 2 dimensional parameter space so we can do the calculation there and plot the result.

In this case, since we are in two dimensions we can use the fft algorithm that is sensibly faster.

In [17]:
param_names = ['delta_omegam', 'delta_sigma8']
shift_P, shift_low, shift_hi = mcmc_tension.kde_parameter_shift_2D_fft(diff_chain, param_names=param_names, feedback=0)

# print the results:
print(f'Shift probability considering all parameters:\n',
      f'   P = {shift_P:.5f} +{shift_hi-shift_P:.5f} -{shift_P-shift_low:.5f}')
# turn the result to effective number of sigmas:
print(f'    n_sigma = {utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'+{utilities.from_confidence_to_sigma(shift_hi)-utilities.from_confidence_to_sigma(shift_P):.3f}',
      f'-{utilities.from_confidence_to_sigma(shift_P)-utilities.from_confidence_to_sigma(shift_low):.3f}')

# triangle plot with the 2D shift probability:
g = plots.get_single_plotter()
diff_chain.updateSettings({'contours': [0.68, 0.95, shift_P]})
g.settings.num_plot_contours = 3
g.triangle_plot(diff_chain, param_names, filled=True, markers={name:0. for name in param_names});
plt.show();

Shift probability considering all parameters:
    P = 0.99769 +0.00004 -0.00004
    n_sigma = 3.047 +0.005 -0.005

<Figure size 600x450 with 0 Axes>
No description has been provided for this image