<jupyter_start><jupyter_text>MadMiner particle physics tutorial Part 4a: Limit setting<NAME>, <NAME>, <NAME>, and <NAME> 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<jupyter_code>from __future__ import absolute_import, division, print_function, unicode_literals

import six
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)<jupyter_output><empty_output><jupyter_text>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`:<jupyter_code>limits = AsymptoticLimits('data/delphes_data_shuffled.h5')<jupyter_output>13:54 madminer.analysis    INFO    Loading data from data/delphes_data_shuffled.h5
13:54 madminer.analysis    INFO    Found 1 parameters
13:54 madminer.analysis    INFO    Did not find nuisance parameters
13:54 madminer.analysis    INFO    Found 3 benchmarks, of which 3 physical
13:54 madminer.analysis    INFO    Found 5 observables
13:54 madminer.analysis    INFO    Found 1661817 events
13:54 madminer.analysis    INFO      567199 signal events sampled from benchmark no-higgs
13:54 madminer.analysis    INFO      533072 signal events sampled from benchmark sm
13:54 madminer.analysis    INFO      561546 signal events sampled from benchmark 5sq-higgs
13:54 madminer.analysis    INFO    Found morphing setup with 3 components
13:54 madminer.analysis    INFO    Did not find nuisance morphing setup<jupyter_text>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.<jupyter_code>grid_ranges = [(0., 5.)]
grid_resolutions = [25]<jupyter_output><empty_output><jupyter_text>What luminosity (in inverse pb) are we talking about?<jupyter_code>lumi = 100.
p_values = {}
mle = {}<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>theta_grid, p_values_expected_histo, best_fit_expected_histo, _, _, (histos, observed, observed_weights) = limits.expected_limits(
    mode="histo",
    hist_vars=["pt_j1"],
    #hist_vars=["delta_eta_jj"],
    theta_true=[5.],
    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<jupyter_output>14:19 madminer.limits      INFO    Generating Asimov data
14:19 madminer.limits      INFO    Expected events: 2344.564619230159
14:19 madminer.limits      INFO    Setting up standard summary statistics
14:19 madminer.limits      INFO    Creating histograms of 1 summary statistics. Using (25,) bins each, or 25 in total.
14:19 madminer.limits      INFO    Determining fixed adaptive histogram binning for all points on grid
14:19 madminer.limits      INFO    Calculating kinematic log likelihood with histograms
14:19 madminer.limits      WARNING Removing %s inf / nan results from calculation
14:19 madminer.limits      INFO    Calculating p-values<jupyter_text>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:<jupyter_code>len(histos)
#indices = [12 + i * 25 for i in [6,9,12,15,18]]
#indices = [6,9,12,15,18]
indices = [0,5,10,15,20]

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.,400.),
    #xrange=(0.,5.),
)

plt.show()<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>theta_grid, p_values_expected_alices, best_fit_expected_alices, _, _, _ = limits.expected_limits(
    mode="ml",
    model_file='models/alices',
    theta_true=[1.],
    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<jupyter_output>14:06 madminer.limits      INFO    Generating Asimov data
14:06 madminer.limits      INFO    Expected events: 24.57890075732583
14:06 madminer.limits      INFO    Loading kinematic likelihood ratio estimator
14:06 madminer.ml          INFO    Loading model from models/alices
14:06 madminer.ml          WARNING Parameter scaling information not found in models/alices
14:06 madminer.limits      INFO    Calculating kinematic log likelihood ratio with estimator
14:06 madminer.ml          INFO    Loading evaluation data
14:06 madminer.ml          INFO    Starting ratio evaluation for 2655500 x-theta combinations
14:07 madminer.ml          INFO    Evaluation done
14:07 madminer.limits      INFO    Calculating p-values<jupyter_text>4. Expected limits based on score estimators To get p-values from a SALLY estimator, we have to use histograms of the estimated score:<jupyter_code>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=[1.],
    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<jupyter_output>14:07 madminer.limits      INFO    Generating Asimov data
14:07 madminer.limits      INFO    Expected events: 24.57890075732583
14:07 madminer.limits      INFO    Loading score estimator and setting all components up as summary statistics
14:07 madminer.ml          INFO    Loading model from models/sally
14:07 madminer.ml          INFO    Starting score evaluation
14:07 madminer.limits      INFO    Creating histograms of 1 summary statistics. Using (25,) bins each, or 25 in total.
14:07 madminer.limits      INFO    Determining fixed adaptive histogram binning for all points on grid
14:07 madminer.ml          INFO    Starting score evaluation
14:07 madminer.ml          INFO    Starting score evaluation
14:07 madminer.limits      INFO    Calculating kinematic log likelihood with histograms
14:07 madminer.limits      WARNING Removing %s inf / nan results from calculation
14:07 madminer.limits      INFO    Calculating p-values<jupyter_text>Let's have a look at the underlying 2D histograms:<jupyter_code>#indices = [12 + i * 25 for i in [0,6,12,18,24]]
indices = [0,5,10,15,20]

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=(0.,.5),
    yrange=(-3.,3.),
    log=True,
    zrange=(1.e-3,1.),
    markersize=10.
)<jupyter_output>/usr/local/lib/python2.7/dist-packages/matplotlib/axes/_base.py:3477: UserWarning: Attempted to set non-positive ylimits for log-scale axis; invalid limits will be ignored.
  'Attempted to set non-positive ylimits for log-scale axis; '<jupyter_text>6. Toy signal In addition to these expected limits (based on the SM), let us inject a mock signal. We first generate the data:<jupyter_code>sampler = SampleAugmenter('data/delphes_data_shuffled.h5')
x_observed, _, _ = sampler.sample_test(
    theta=sampling.morphing_point([4.]),
    n_samples=1000,
    folder=None,
    filename=None,
)
obs_n = 79
_, p_values_observed, best_fit_observed, _, _, _ = limits.observed_limits(
    n_observed=obs_n,
    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 signal"] = p_values_observed
mle["ALICES signal"] = best_fit_observed
theta_grid, p_values_expected_alices, best_fit_expected_alices, _, _, _ = limits.expected_limits(
    mode="ml",
    model_file='models/alices',
    theta_true=[4.],
    grid_ranges=grid_ranges,
    grid_resolutions=grid_resolutions,
    luminosity=lumi,
    include_xsec=False,
)

p_values["ALICES signal"] = p_values_expected_alices
mle["ALICES signal"] = best_fit_expected_alices
theta_grid, p_values_expected_sally, best_fit_expected_sally, _, _, _ = limits.observed_limits(
    mode="sally",
    model_file='models/sally',
    n_observed=obs_n,
    x_observed=x_observed,
    grid_ranges=grid_ranges,
    grid_resolutions=grid_resolutions,
    luminosity=lumi,
    include_xsec=False,
)

p_values["SALLY signal"] = p_values_expected_sally
mle["SALLY signal"] = best_fit_expected_sally
theta_grid, p_values_expected_histo, best_fit_expected_histo, _, _, _ = limits.observed_limits(
    mode="histo",
    hist_vars=["pt_j1"],
    n_observed=obs_n,
    x_observed=x_observed,
    grid_ranges=grid_ranges,
    grid_resolutions=grid_resolutions,
    luminosity=lumi,
    include_xsec=False,
)

p_values["Histogram signal"] = p_values_expected_histo
mle["Histogram signal"] = best_fit_expected_histo
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=[4.],
    grid_ranges=grid_ranges,
    grid_resolutions=grid_resolutions,
    luminosity=lumi,
    include_xsec=False,
    return_asimov=True
)

p_values["Histogram signal"] = p_values_expected_histo
mle["Histogram signal"] = best_fit_expected_histo<jupyter_output>14:08 madminer.limits      INFO    Generating Asimov data
14:08 madminer.limits      INFO    Expected events: 1413.9806174826467
14:08 madminer.limits      INFO    Setting up standard summary statistics
14:08 madminer.limits      INFO    Creating histograms of 1 summary statistics. Using (25,) bins each, or 25 in total.
14:08 madminer.limits      INFO    Determining fixed adaptive histogram binning for all points on grid
14:08 madminer.limits      INFO    Calculating kinematic log likelihood with histograms
14:08 madminer.limits      WARNING Removing %s inf / nan results from calculation
14:08 madminer.limits      INFO    Calculating p-values<jupyter_text>Plot Let's plot the results:<jupyter_code>show = "ALICES"
cmin, cmax = 1.e-3, 1.

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])

for i, (label, p_value) in enumerate(six.iteritems(p_values)):
    if (label.find("signal") != -1):
        continue
    #plt.scatter(centers**2, p_value, label=label)
    plt.scatter(centers, p_value, label=label)
    print (theta_grid[mle[label]], label)
#     plt.scatter(
#         theta_grid[mle[label]], theta_grid[mle[label]][1],
#         s=80., color='C{}'.format(i), marker='*',
#         label=label
#     )

plt.legend()

plt.xlabel(r'$\theta_0$')
#plt.xlabel(r'$\mu$')
plt.ylabel(r'p_value')
#plt.xlim(0,5)

plt.tight_layout()
plt.show()


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])

for i, (label, p_value) in enumerate(six.iteritems(p_values)):
    if (label.find("signal") == -1):
        continue
    plt.scatter(centers, p_value, label=label)
    print (theta_grid[mle[label]], label)
#     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'p_value')

plt.tight_layout()
plt.show()<jupyter_output>[0.83333333] ALICES signal
[0.625] Histogram signal
{'fraction_non_alphanumeric': 0.0985414701400479, 'fraction_numerical': 0.037587983455482185, 'mean_word_length': 4.602439024390244, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '46261', 'license': 'MIT', 'n_tokens_mistral': 4938, 'n_tokens_neox': 4511, 'n_words': 1189, 'path': 'examples/tutorial_h4l/defunct_4a_limits_JB.ipynb', 'repo_name': 'aghoshpub/LikelihoodFreeInterference', 'size': '13806'}