MadMiner physics tutorial (part 4A)

Johann Brehmer, Felix Kling, Irina Espejo, and Kyle Cranmer 2018-2019

In part 4A of this tutorial we will use the networks trained in step 3A and 3B to calculate the expected limits on our theory parameters.

0. Preparations

import logging
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline

from madminer.limits import AsymptoticLimits
from madminer.sampling import SampleAugmenter
from madminer import sampling
from madminer.plotting import plot_histograms
# MadMiner output
logging.basicConfig(
    format='%(asctime)-5.5s %(name)-20.20s %(levelname)-7.7s %(message)s',
    datefmt='%H:%M',
    level=logging.INFO
)

# Output of all other modules (e.g. matplotlib)
for key in logging.Logger.manager.loggerDict:
    if "madminer" not in key:
        # print("Deactivating logging output for", key)
        logging.getLogger(key).setLevel(logging.WARNING)

1. Preparations

In the end, what we care about are not plots of the log likelihood ratio, but limits on parameters. But at least under some asymptotic assumptions, these are directly related. MadMiner makes it easy to calculate p-values in the asymptotic limit with the AsymptoticLimits class in the madminer.limits:

limits = AsymptoticLimits('data/lhe_data_shuffled.h5')
# limits = AsymptoticLimits('data/delphes_data_shuffled.h5')
20:13 madminer.analysis    INFO    Loading data from data/lhe_data_shuffled.h5
20:13 madminer.analysis    INFO    Found 2 parameters
20:13 madminer.analysis    INFO    Did not find nuisance parameters
20:13 madminer.analysis    INFO    Found 6 benchmarks, of which 6 physical
20:13 madminer.analysis    INFO    Found 3 observables
20:13 madminer.analysis    INFO    Found 89991 events
20:13 madminer.analysis    INFO      49991 signal events sampled from benchmark sm
20:13 madminer.analysis    INFO      10000 signal events sampled from benchmark w
20:13 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_w
20:13 madminer.analysis    INFO      10000 signal events sampled from benchmark ww
20:13 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_ww
20:13 madminer.analysis    INFO    Found morphing setup with 6 components
20:13 madminer.analysis    INFO    Did not find nuisance morphing setup

This class provids two high-level functions:

  • AsymptoticLimits.observed_limits() lets us calculate p-values on a parameter grid for some observed events, and

  • AsymptoticLimits.expected_limits() lets us calculate expected p-values on a parameter grid based on all data in the MadMiner file.

First we have to define the parameter grid on which we evaluate the p-values.

grid_ranges = [(-1, 1.), (-1, 1.)]
grid_resolutions = [25, 25]

What luminosity (in inverse pb) are we talking about?

lumi = 300.*1000. 

2. Expected limits based on histogram

First, as a baseline, let us calculate the expected limits based on a simple jet pT histogram. Right now, there are not a lot of option for this function; MadMiner even calculates the binning automatically. (We will add more functionality!)

The keyword include_xsec determines whether we include information from the total rate or just use the shapes. Since we don’t model backgrounds and systematics in this tutorial, the rate information is unrealistically large, so we leave it out here.

## create dictionaries to store results
p_values = {}
mle = {}
try:
    theta_grid, p_values_expected_histo, best_fit_expected_histo, _, _, (histos, observed, observed_weights) = limits.expected_limits(
        mode="histo",
        hist_vars=["pt_j1"],
        theta_true=[0.,0.],
        grid_ranges=grid_ranges,
        grid_resolutions=grid_resolutions,
        luminosity=lumi,
        include_xsec=False,
        return_asimov=True,
    )

    p_values["Histogram"] = p_values_expected_histo
    mle["Histogram"] = best_fit_expected_histo
except: 
    print("no Histogram")
20:13 madminer.limits      INFO    Generating Asimov data
20:13 madminer.limits      INFO    Expected events: 296.8608606876634
20:13 madminer.limits      INFO    Calculating p-values for 296.8608606876634 expected events in mode histo without rate information
20:13 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:13 madminer.limits      INFO    Setting up standard summary statistics
20:13 madminer.limits      INFO    Creating histograms of 1 summary statistics. Using (25,) bins each, or 25 in total.
20:13 madminer.limits      INFO    Determining fixed adaptive histogram binning for all points on grid
20:13 madminer.analysis    WARNING Requested 100000 events, but only 71993 available
20:13 madminer.analysis    WARNING Requested 100000 events, but only 71993 available
20:13 madminer.limits      INFO    Calculating kinematic log likelihood with histograms
20:13 madminer.limits      INFO    Calculating p-values

With mode="rate", we could calculate limits based on only the rate – but again, since the rate is extremely powerful when backgrounds and systematics are not taken into account, we don’t do that in this tutorial.

Let’s visualize the likelihood estimated with these histograms:

indices = [12 + i * 25 for i in [6,9,12,15,18]]

fig = plot_histograms(
    histos=[histos[i] for i in indices],
    observed=[observed[i] for i in indices],
    observed_weights=observed_weights,
    histo_labels=[r"$\theta_0 = {:.2f}$".format(theta_grid[i,0]) for i in indices],
    xlabel="Jet $p_T$",
    xrange=(0.,500.),
)

plt.show()
../../../_images/4a_limits_18_0.png

3. Expected limits based on ratio estimators

Next, mode="ml" allows us to calculate limits based on any ParamterizedRatioEstimator instance like the ALICES estimator trained above:

try:
    theta_grid, p_values_expected_alices, best_fit_expected_alices, _, _, _ = limits.expected_limits(
        mode="ml",
        model_file='models/alices',
        theta_true=[0.,0.],
        grid_ranges=grid_ranges,
        grid_resolutions=grid_resolutions,
        luminosity=lumi,
        include_xsec=False,
    )

    p_values["ALICES"] = p_values_expected_alices
    mle["ALICES"] = best_fit_expected_alices
except:
    print("no ALICES")
20:13 madminer.limits      INFO    Generating Asimov data
20:13 madminer.limits      INFO    Expected events: 296.8608606876634
20:13 madminer.limits      INFO    Calculating p-values for 296.8608606876634 expected events in mode ml without rate information
20:13 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:13 madminer.limits      INFO    Loading kinematic likelihood ratio estimator
20:13 madminer.ml          INFO    Loading model from models/alices
20:13 madminer.limits      INFO    Calculating kinematic log likelihood ratio with estimator
20:13 madminer.ml          INFO    Loading evaluation data
20:13 madminer.ml          INFO    Starting ratio evaluation for 6310625 x-theta combinations
20:16 madminer.ml          INFO    Evaluation done
20:16 madminer.limits      INFO    Calculating p-values

and the ALICES estimator which only used the observable pt_j1

try:
    theta_grid, p_values_expected_alices_pt, best_fit_expected_alices_pt, _, _, _ = limits.expected_limits(
        mode="ml",
        model_file='models/alices_pt',
        theta_true=[0.,0.],
        grid_ranges=grid_ranges,
        grid_resolutions=grid_resolutions,
        luminosity=lumi,
        include_xsec=False,
    )

    p_values["ALICES_pt"] = p_values_expected_alices_pt
    mle["ALICES_pt"] = best_fit_expected_alices_pt
except:
    print("no ALICES_pt")
20:16 madminer.limits      INFO    Generating Asimov data
20:16 madminer.limits      INFO    Expected events: 296.8608606876634
20:16 madminer.limits      INFO    Calculating p-values for 296.8608606876634 expected events in mode ml without rate information
20:16 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:16 madminer.limits      INFO    Loading kinematic likelihood ratio estimator
20:16 madminer.ml          INFO    Loading model from models/alices_pt
20:16 madminer.limits      INFO    Calculating kinematic log likelihood ratio with estimator
20:16 madminer.ml          INFO    Loading evaluation data
20:16 madminer.ml          INFO    Starting ratio evaluation for 6310625 x-theta combinations
20:18 madminer.ml          INFO    Evaluation done
20:18 madminer.limits      INFO    Calculating p-values

4. Expected limits based on score estimators

To get p-values from a SALLY estimator, we have to use histograms of the estimated score:

try:
    theta_grid, p_values_expected_sally, best_fit_expected_sally, _, _, (histos, observed, observed_weights) = limits.expected_limits(
        mode="sally",
        model_file='models/sally',
        theta_true=[0.,0.],
        grid_ranges=grid_ranges,
        grid_resolutions=grid_resolutions,
        luminosity=lumi,
        include_xsec=False,
        return_asimov=True,
    )
    p_values["SALLY"] = p_values_expected_sally
    mle["SALLY"] = best_fit_expected_sally
except: 
    print("no SALLY")
20:18 madminer.limits      INFO    Generating Asimov data
20:18 madminer.limits      INFO    Expected events: 296.8608606876634
20:18 madminer.limits      INFO    Calculating p-values for 296.8608606876634 expected events in mode sally without rate information
20:18 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:18 madminer.limits      INFO    Loading score estimator and setting all components up as summary statistics
20:18 madminer.ml          INFO    Loading model from models/sally
20:18 madminer.ml          INFO    Starting score evaluation
20:18 madminer.limits      INFO    Creating histograms of 2 summary statistics. Using (8, 8) bins each, or 64 in total.
20:18 madminer.limits      INFO    Determining fixed adaptive histogram binning for all points on grid
20:18 madminer.analysis    WARNING Requested 100000 events, but only 71993 available
20:18 madminer.ml          INFO    Starting score evaluation
20:18 madminer.analysis    WARNING Requested 100000 events, but only 71993 available
20:18 madminer.ml          INFO    Starting score evaluation
20:19 madminer.limits      INFO    Calculating kinematic log likelihood with histograms
20:19 madminer.limits      INFO    Calculating p-values

Let’s have a look at the underlying 2D histograms:

indices = [12 + i * 25 for i in [0,6,12,18,24]]

fig = plot_histograms(
    histos=[histos[i] for i in indices],
    observed=observed[0,:100,:],
    observed_weights=observed_weights[:100],
    histo_labels=[r"$\theta_0 = {:.2f}$".format(theta_grid[i,0]) for i in indices],
    xlabel=r'$\hat{t}_0(x)$',
    ylabel=r'$\hat{t}_1(x)$',
    xrange=(-3.,.5),
    yrange=(-3.,3.),
    log=True,
    zrange=(1.e-3,1.),
    markersize=10.
)
    
../../../_images/4a_limits_28_0.png

5. Expected limits based on likelihood estimators

try:
    theta_grid, p_values_expected_scandal, best_fit_expected_scandal, _, _, _ = limits.expected_limits(
        mode="ml",
        model_file='models/scandal',
        theta_true=[0.,0.],
        grid_ranges=grid_ranges,
        grid_resolutions=grid_resolutions,
        luminosity=lumi,
        include_xsec=False,
    )

    p_values["SCANDAL"] = p_values_expected_scandal
    mle["SCANDAL"] = best_fit_expected_scandal
except:
    print("no SCANDAL")
20:19 madminer.limits      INFO    Generating Asimov data
20:19 madminer.limits      INFO    Expected events: 296.8608606876634
20:19 madminer.limits      INFO    Calculating p-values for 296.8608606876634 expected events in mode ml without rate information
20:19 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:19 madminer.limits      INFO    Loading kinematic likelihood ratio estimator
no SCANDAL

6. Toy signal

In addition to these expected limits (based on the SM), let us inject a mock signal. We first generate the data:

sampler = SampleAugmenter('data/lhe_data_shuffled.h5')
x_observed, _, _ = sampler.sample_test(
    theta=sampling.morphing_point([0.5,0.5]),
    n_samples=100,
    folder=None,
    filename=None,
)
20:19 madminer.analysis    INFO    Loading data from data/lhe_data_shuffled.h5
20:19 madminer.analysis    INFO    Found 2 parameters
20:19 madminer.analysis    INFO    Did not find nuisance parameters
20:19 madminer.analysis    INFO    Found 6 benchmarks, of which 6 physical
20:19 madminer.analysis    INFO    Found 3 observables
20:19 madminer.analysis    INFO    Found 89991 events
20:19 madminer.analysis    INFO      49991 signal events sampled from benchmark sm
20:19 madminer.analysis    INFO      10000 signal events sampled from benchmark w
20:19 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_w
20:19 madminer.analysis    INFO      10000 signal events sampled from benchmark ww
20:19 madminer.analysis    INFO      10000 signal events sampled from benchmark neg_ww
20:19 madminer.analysis    INFO    Found morphing setup with 6 components
20:19 madminer.analysis    INFO    Did not find nuisance morphing setup
20:19 madminer.sampling    INFO    Extracting evaluation sample. Sampling according to [0.5 0.5]
20:19 madminer.sampling    INFO    Starting sampling serially
20:19 madminer.sampling    INFO    Sampling from parameter point 1 / 1
20:19 madminer.sampling    INFO    Effective number of samples: mean 195.22279009818115, with individual thetas ranging from 195.22279009818112 to 195.22279009818112
_, p_values_observed, best_fit_observed, _, _, _ = limits.observed_limits(
    x_observed=x_observed,
    mode="ml",
    model_file='models/alices',
    grid_ranges=grid_ranges,
    grid_resolutions=grid_resolutions,
    luminosity=lumi,
    include_xsec=False,
)

p_values["ALICES (observed)"] = p_values_observed
mle["ALICES (observed)"] = best_fit_observed
20:19 madminer.limits      INFO    Calculating p-values for 100 expected events in mode ml without rate information
20:19 madminer.limits      INFO    Evaluating likelihood on a regular grid with 625 parameter points
20:19 madminer.limits      INFO    Loading kinematic likelihood ratio estimator
20:19 madminer.ml          INFO    Loading model from models/alices
20:19 madminer.limits      INFO    Calculating kinematic log likelihood ratio with estimator
20:19 madminer.ml          INFO    Loading evaluation data
20:19 madminer.ml          INFO    Starting ratio evaluation for 62500 x-theta combinations
20:19 madminer.ml          INFO    Evaluation done
20:19 madminer.limits      INFO    Calculating p-values

7. Plot

Let’s plot the results:

show = "ALICES"

bin_size = (grid_ranges[0][1] - grid_ranges[0][0])/(grid_resolutions[0] - 1)
edges = np.linspace(grid_ranges[0][0] - bin_size/2, grid_ranges[0][1] + bin_size/2, grid_resolutions[0] + 1)
centers = np.linspace(grid_ranges[0][0], grid_ranges[0][1], grid_resolutions[0])

fig = plt.figure(figsize=(6,5))
ax = plt.gca()

cmin, cmax = 1.e-2, 1.

pcm = ax.pcolormesh(
    edges, edges, p_values[show].reshape((grid_resolutions[0], grid_resolutions[1])).T,
    norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
    cmap='Greys_r'
)
cbar = fig.colorbar(pcm, ax=ax, extend='both')

for i, (label, p_value) in enumerate(p_values.items()):
    plt.contour(
        centers, centers, p_value.reshape((grid_resolutions[0], grid_resolutions[1])).T,
        levels=[0.32],
        linestyles='-', colors='C{}'.format(i)
    )
    plt.scatter(
        theta_grid[mle[label]][0], theta_grid[mle[label]][1],
        s=80., color='C{}'.format(i), marker='*',
        label=label
    )

plt.legend()

plt.xlabel(r'$\theta_0$')
plt.ylabel(r'$\theta_1$')
cbar.set_label('Expected p-value ({})'.format(show))

plt.tight_layout()
plt.show()
../../../_images/4a_limits_37_0.png

Save results

Let’s save the results to later compare them to results obtained in other notebooks.

np.save("limits/limits.npy", [p_values,mle])