<jupyter_start><jupyter_text>Making a Combined Fit with the Asymptotic Formulas In this exercise we will go beyond estimating epected significances with the **approximate median significance** calculated for each bin assuming Poissonian counts, and then adding them in quadrature.For a complete analysis, this approach that sees all bins as uncorrelated is too simplistic and will become invalid as soon as you introduce things like systematic uncertainties which are correlated between the bins.Therefore, we need to set up a likelihood-based fitting framework which is a bit more flexible. If you are curious how the CMS collaboration does this, keep in mind that it has it's own framework called "[Higgs combine](https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit)". We won't use it because it is a big "black box" and does not have much pedagogical value.The only black box we will use for our fits is the ancient but proven [MINUIT](https://en.wikipedia.org/wiki/MINUIT) minimization package (which is also used by [ROOT](https://root.cern.ch/)). I particular, we use it's python wrapper [iminuit](https://github.com/scikit-hep/iminuit). Of course there are other minimation package in python that might be even more wide-spread (the most famous one is [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html), which wraps around the Fortran code [MINPACK](https://en.wikipedia.org/wiki/MINPACK)).In the `minuit_minimize.py` file, I implemented a function that uses *iminuit* to minimize a function that takes a dictionary of parameters. Let's import it, and of couse *pandas* and *numpy* too.<jupyter_code>import numpy as np
import pandas as pd

from geeksw.optimize.minuit_minimize import minuit_minimize<jupyter_output><empty_output><jupyter_text>If you want to see how that function can be used, run `?minuit_minimize` or `??minuit_minimize` to inspect the documentation string or source code.Next, we summarize our results in a data frame that we will then pass to some fitting function.In this example, we will use some values from an outdated version of the CMS analysis, which also includes more final states than just the four leptons.The rows will correspond to the category bins that we have in the analysis. In this example, the bins are the  different lepton multiplicities in the final state (`SSplusJets` means two same sign leptons plus jets and is used in the *WWW* analysis).The columns are the different components that we have in the simulation. At this point it doesn't matter how the backgrounds are composed, and we just have one background component `B` in addition to four on-shell *VVV* and four *VH* to *VVV* components.<jupyter_code>datacard = pd.DataFrame(dict(WWW=[27.90, 16.07, 0, 0, 0],
                             WWZ=[0.91, 0.91, 8.57, 0, 0],
                             WZZ=[0.18, 0.01, 0.34, 1.08, 0.01],
                             ZZZ=[0, 0, 0.07, 0.25, 0.41],
                             VH_WWW=[17.62, 8.47, 0, 0, 0],
                             VH_WWZ=[0, 0,  5.6+0.09+0.17, 0, 0],
                             VH_WZZ=[0.0, 0, 0, 0, 0],
                             VH_ZZZ=[0.0, 0, 0, 0, 0],
                             B=[447.21, 94.99, 15.76, 0.55, 0.06]),
                        index=["SSplusJets", "3leptons", "4leptons", "5leptons", "6leptons"])
datacard.index.name = "bin"<jupyter_output><empty_output><jupyter_text>Let's check out the data frame, which we calll "datacard", because that's how we call these kind of summary tables that we put into fitting tools like *Higgs combine* inside CMS:<jupyter_code>datacard<jupyter_output><empty_output><jupyter_text>Now it will get interesting. We want to use the recommended statistical analysis for LHC experiments described in [this paper](https://arxiv.org/pdf/1503.07622.pdf), in particular Chapter 5 about frequentist statistical procedures.In a discovery setting where we don't know a priori the signal strength $\mu$ (the ratio of observed over expected signal), we compute the likelihood ratio to test the background hypothesis:$$ \lambda(0) = \frac{L(0, \hat{\hat{\theta}}(0))}{L(\hat{\mu}, \hat{\theta})}. $$This corresponds to Equation 49 in the paper with the expected signal $\mu = 0$. Next, we compute the discovery test statistic $q_0$ like in Equation 51:$$q_0 = -2\log{\lambda(0)}.$$Finally, we use the asymptotic formula from Equation 72 to get the significance:$$Z = \sqrt{q_0}.$$In the following, I show you my implementation of the LHC discovery test. It starts with some imports and a little helper function which is not important to understand now:<jupyter_code>import inspect

import scipy.special

def get_nuisance_parameter_name(func):
    args = inspect.getfullargspec(func).args
    assert(len(args) == 2)
    assert(args[1] == "datacard")
    return args[0]<jupyter_output><empty_output><jupyter_text>Next, I have a function that builds **negative log-likelihood functions** which you have to pass:* a datacard* the definition of signal strengths we want to measure* a list of nuisance parameter functions and penalties (we will cover this later)* the number of observed events in data if you want to compute the observed significance (if this is `None` you will get the expected significance)<jupyter_code>def make_nll_func(datacard, signal_strengths, nuisances=[], observations=None):
        
    if observations is None:
        observations = datacard.sum(axis=1)
    
    # This is the negative logarithm of the Possonian probabiliy density to observe n events
    # with an expectation s + b
    def nll_term(s, b, n):
        return (s+b) - n * np.log(s+b)
    
    signal_columns = "+".join(signal_strengths.values()).split("+")
    background_columns = [c for c in datacard.columns if not c in signal_columns]

    nuisance_param_names = [get_nuisance_parameter_name(nuis) for nuis in nuisances]
    
    def nll_func(**params):
                
        df = datacard.copy(deep=True)
        
        nll = 0.0

        for nuisance, param_name in zip(nuisances, nuisance_param_names):
            df, constraint = nuisance(params[param_name], df)
            nll += constraint
        
        s = pd.DataFrame({key : params[key] * \
                          df.eval(signal_strengths[key]) for key in signal_strengths}).sum(axis=1)
        b = df[background_columns].sum(axis=1)
    
        return nll + nll_term(s, b, observations).sum()
    
    return nll_func<jupyter_output><empty_output><jupyter_text>Having the full **negative log-likelihood** function for you analysis is very valuable. Next, we will use it to compute discovery significances in this `asymptotic_discovery_fit` function. It takes the same arguments as `make_nll_func` but additionally **iminuit** to find optimal parameter values and execute the LHC discovery statistic prescription.In the resulting tuple, you will have one data frame which shows you the fitted signal strengths and discovery significances, and one data frame which summarizes the nuisance parameter fit result (not important for now).<jupyter_code>def asymptotic_discovery_fit(datacard, signal_strengths=None, nuisances=[], observations=None):

    nll_func = make_nll_func(datacard, signal_strengths, nuisances, observations)
 
    nuisance_param_names = [get_nuisance_parameter_name(nuis) for nuis in nuisances]
    signal_strength_param_names = list(signal_strengths.keys())
    param_names = signal_strength_param_names + nuisance_param_names
    n_params = len(param_names)
    starting_params = [1.1] * n_params
    starting_errors = [0.1] * n_params
    
    values, errors, _ = minuit_minimize(nll_func, param_names, starting_params, starting_errors)
    
    df = pd.DataFrame(dict(value=list(values.values()),
                           error=list(errors.values())), index=param_names)
    df = df.eval("precision=error/value")
    
    def set_elem_to_zero(a,  name):
        b = dict(**a)
        b[name] = 0.0
        return b
    
    def two_dnnl(name):
        fixed_params = {name : 0.0}
        if name in fixed_params and len(param_names) == 1:
            null_nll = nll_func(**fixed_params)
        else:
            null_nll = minuit_minimize(nll_func, param_names, starting_params, starting_errors,
                                       fixed_params=fixed_params).f_val
        ll_0 = -null_nll
        ll_best_fit = -nll_func(**values)
        return -2 * (ll_0 - ll_best_fit)

    df_mu = df.loc[signal_strength_param_names]
    df_nuisance = df.loc[nuisance_param_names]
    
    df_mu["significance"] = np.sqrt(np.array(list(map(two_dnnl, signal_strength_param_names))))
    df_mu["p_value"] = 0.5 * (1.+scipy.special.erf(-df_mu["significance"]/np.sqrt(2.)))
    
    return df_mu, df_nuisance<jupyter_output><empty_output><jupyter_text>I will show you how to use this function. The **signal strength** definitions are just a dictionary of data frame queries, and each entry is a separate signal for the analysis. Therefore, we can very quickly try out different signal definitions.For example, taking the sum of all triboson final states as one single signal:<jupyter_code>signal_strengths = dict(mu_all="WWW+WWZ+WZZ+ZZZ+VH_WWW+VH_WWZ+VH_WZZ+VH_ZZZ")<jupyter_output><empty_output><jupyter_text>Or consider only on-shell production and have separate signal strengths for different final states:<jupyter_code>signal_strengths = dict(mu_WWW="WWW", mu_WWZ="WWZ", mu_WZZ="WZZ", mu_ZZZ="ZZZ")<jupyter_output><empty_output><jupyter_text>Pass the datacard and the signal strength definition to the function and save the results:<jupyter_code>df_mu, df_nuisance = asymptotic_discovery_fit(datacard, signal_strengths)<jupyter_output><empty_output><jupyter_text>Let's look at the results (`df_nuisance` is empty if you didn't pass nuisance parameters):<jupyter_code>df_mu<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09057639524245197, 'fraction_numerical': 0.015655179424621326, 'mean_word_length': 3.8630746416213544, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '426842', 'license': 'MIT', 'n_tokens_mistral': 3113, 'n_tokens_neox': 2897, 'n_words': 1024, 'path': 'examples/stats/05_AsymptoticSignificanceFit.ipynb', 'repo_name': 'guitargeek/geeksw', 'size': '9837'}